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Abstract 


Concurrent  engineering  approaches  for  the  disciplines  of  computational  fluid  dynamics 
(CFD)  and  electromagnetics  (CEM)  are  necessary  for  designing  future  high-performance,  low- 
observable  aircraft.  A  characteristic-based  finite-volume  time-domain  (FVTD)  computational 
algorithm  developed  for  CFD  and  herein  applied  to  CEM  is  implemented  to  analyze  the  radar 
cross  section  (RCS)  of  two  three-dimensional  objects,  the  ogive  and  cone-sphere,  by  utilizing  a 
scattered-field  formulation  of  the  time-dependent  Maxwell  equations.  The  FVTD  formulation 
uses  a  van  Leer’ s  kappa  scheme  for  the  flux  evaluation  and  a  Runge-Kutta  multi-stage  scheme 
for  the  time  integration.  The  RCS  results  are  obtained  from  the  electromagnetic  fields 
subsequent  to  a  Fourier  transform  and  a  near-to-far  field  transformation. 

Developmental  work  for  the  thesis  focused  on  modifying  the  original  code  to  analyze 
scattering  and  obtain  RCS  data  for  closed-surface  perfect  electric  conductor  (PEC)  3-D  objects 
using  either  a  Gaussian  pulse  or  sinusoid  incident  wave.  Specification  of  the  direction  and 
polarization  of  the  incident  wave  provides  monostatic  and  bistatic  results  for  various 
polarizations.  A  RCS  convergence  check,  used  with  a  sinusoid,  ends  the  simulation  after  the 
transients  diminish  and  the  bistatic  RCS  values  are  within  0. 1  dB  of  the  RCS  values  calculated 
during  the  previous  period.  A  threshold  check,  used  with  a  Gaussian  pulse,  ends  the  simulation 
once  the  amplitude  of  the  scattered  field  is  140  dB  below  the  peak  of  the  Gaussian  pulse.  A 
bistatic-to-monostatic  RCS  approximation  saves  computer  run  time  by  a  factor  of  nearly  40. 

The  FVTD  code  and  algorithm  are  validated  for  electromagnetic  scattering  problems  by 
comparing  FVTD  code  RCS  results  to  data  obtained  from  the  Moment  Method  (MoM)  code, 
CICERO,  and  empirical  RCS  data  published  by  the  Electromagnetic  Code  Consortium  (EMCC). 
The  FVTD  RCS  results  for  the  ogive  and  cone-sphere  are  within  3.0  dB  of  the  bistatic  MoM 
results  and  3.1  dB  of  the  monostatic  empirical  RCS  data.  Accurate  FVTD  computations  of 
diffraction,  traveling  waves,  and  creeping  waves  require  a  surface  grid  point  density  of  15-30 
cells/?i  and  a  PEC  boundary  condition  grid  density  of  200-400  cells/A.,  dependent  on  frequency. 
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APPLICATION  OF  A  FINITE-VOLUME  TIME- 
DOMAIN  MAXWELL  EQUATION  SOLVER  TO 
THREE-DIMENSIONAL  OBJECTS 

1  Introduction 

When  designing  low-observable,  high-performance  military  aircraft,  engineers  must 
consider  both  a  low  radar  cross  section  (RCS)  and  excellent  aerodynamic  performance. 
However,  an  optimum  electromagnetic  design  may  not  be  an  optimum  aerodynamic  design.  The 
efficient  design  of  future  military  aircraft  requires  concurrent  multi-disciplinary  approaches  for 
the  electromagnetic  and  fluid  dynamic  disciplines.  One  computational  design  tool,  the  finite- 
volume  time-domain  (FVTD)  technique,  can  potentially  consider  optimum  electromagnetic  and 
aerodynamic  designs  concurrently,  benefiting  the  design  process  of  low-observable,  high- 
performance  aircraft  [47].  Researchers  have  used  FVTD  successfully  in  the  computational  fluid 
dynamics  (CFD)  discipline  since  the  early  1980’s  [48]  to  analyze  the  airflow  about  an  aircraft, 
and  recently,  Shang,  Shankar,  Blake,  Bishop,  Huh,  Noack,  and  others  have  applied  the  technique 
to  the  computational  electromagnetics  (CEM)  discipline  [5,  6-7,  15,  22,  29-47].  The  application 
of  the  FVTD  technique  to  the  time-domain  Maxwell  equations  of  electromagnetics  is  relatively 
new  compared  to  other  CEM  techniques  but  is  proving  to  be  successful. 

FVTD  requires  application  to  complicated  three-dimensional  (3-D)  benchmark  objects 
for  researchers  to  consider  it  a  feasible  computational  tool  for  electromagnetic  scattering 
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problems.  In  this  thesis  research,  the  comparison  of  FVTD  scattering  results  for  the  three- 
dimensional  objects,  the  ogive  and  cone-sphere,  to  Method  of  Moments  (MoM)  and  empirical 
RCS  data  validates  the  FVTD  technique  in  the  area  of  electromagnetics.  The  MoM  calculations, 
obtained  from  the  code  CICERO,  and  the  empirical  RCS  data,  published  by  the  Electromagnetic 
Code  Consortium  (EMCC),  serve  as  benchmarks  for  the  FVTD  RCS  computations. 

1.1  Background 

FVTD  is  a  logical  choice  for  simulations  in  both  computational  fluid  dynamics  and 
computational  electromagnetics.  Fluid  dynamics  utilizes  FVTD  to  solve  the  Euler  equations  by 
calculating  the  fluxes,  such  as  energy  or  mass,  through  finite-volume  cells  of  a  discretized  space 
[8].  The  Maxwell  equations,  which  model  electromagnetic  phenomena,  are  hyperbolic  like  the 
Euler  equations;  that  is,  both  systems  of  partial  differential  equations  have  real  eigenvalues 
(characteristic  values).  Their  similar  mathematical  form  permits  the  use  of  the  same  algorithm  to 
numerically  solve  each  set  of  equations.  Thus,  hyperbolic  equation  solvers  developed  for  CFD, 
such  as  FVTD,  can  potentially  be  applied  to  Maxwell’s  equations  of  CEM. 

In  addition  to  multi-discipline  applications,  FVTD  has  other  potential  advantages 
including  [18,  51] 

•  the  direct  solution  of  Maxwell’s  two  curl  partial  differential  equations 

•  the  analysis  of  electromagnetic  propagation  through  frequency-dependent,  time- 
dependent,  and  anisotropic  materials 

•  the  ability  to  obtain  multiple  frequency  results  with  a  single  simulation 
Specifically,  the  direct  solution  of  Maxwell’s  time-dependent  equations  is  robust  and  potentially 
as  accurate  as  the  Method  of  Moments  [18,  51].  The  finite-difference  time-domain  (FDTD) 
technique,  which  also  directly  solves  Maxwell’s  partial  differential  equations,  provides  highly 
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accurate  simulations  for  free  space  electromagnetic  scattering  problems  [51].  The  Helmholtz 
equations,  the  time-harmonic  wave  equations  derived  from  the  Maxwell  equations,  are  not  used; 
therefore,  the  propagation  of  electromagnetic  waves  through  frequency-dependent  and  time- 
dependent  materials  can  be  calculated.  FVTD  also  permits  the  analysis  of  anisotropic  materials, 
such  as  radar  absorbing  materials  (RAM)  found  in  filled  honeycomb  structures  [42]. 
Additionally,  the  time-domain  algorithm  allows  the  solution  of  a  problem  for  more  than  one 
frequency.  Broadband  incident  waves,  such  as  Gaussian  pulses,  can  be  used  to  obtain  multiple 
frequency  results. 

The  two  most  prominent  FVTD  researchers  in  the  area  of  CEM  are  Dr.  Vijaya  Shankar 
of  the  Rockwell  International  Science  Center  [42-47]  and  Dr.  Joe  Shang  of  the  Air  Force’s 
Wright  Laboratory  [29-41].  Both  have  written  characteristic-based  FVTD  codes  for  analyzing 
the  electromagnetic  scattering  from  objects.  The  two  codes  differ  in  several  significant  areas: 
mathematical  algorithm,  order  of  accuracy,  and  capability.  Shang  has  implemented  a  spatially 
third-order  accurate  algorithm  and  a  temporal  fourth-order  accurate  scheme  in  the  FVTD  code, 
while  Shankar  has  implemented  a  second-order  accurate  algorithm.  The  higher-order  accurate 
code  can  potentially  reduce  the  required  number  of  grid  points  per  wavelength  permitting  the 
computation  of  the  electromagnetic  scattering  from  electrically  larger  objects  and  reducing  the 
computer  simulation  time  [36,  40].  The  differences  between  the  algorithms  and  the  accuracy  of 
the  codes  are  discussed  in  more  detail  in  Chapter  2. 

The  Rockwell  FVTD  code  currently  has  more  capability  than  the  Wright  Laboratory 
FVTD  code.  The  Rockwell  FVTD  code  has  the  ability  to  analyze  the  scattering  from 
complicated  objects  and  surfaces  such  as  airfoils,  inlets,  and  edges,  while  the  Wright  Lab  FVTD 
code,  prior  to  this  research,  was  limited  to  calculating  the  scattering  from  a  sphere  [39,  41,  44, 
46].  The  Rockwell  code  can  also  calculate  the  electromagnetic  scattering  from  layered  dielectric 
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surfaces  such  as  those  found  on  low-observable  aircraft.  The  Wright  Lab  FVTD  requires 


modification  to  calculate  the  electromagnetic  scattering  from  complicated  3-D  objects  and 
dielectric  surfaces.  After  first  modifying  the  Wright  Lab  FVTD  code  for  this  thesis  research,  an 
analysis  of  the  electromagnetic  scattering  is  performed  for  the  complicated  perfect  electric 
conductor  (PEC)  objects:  ogive  and  cone-sphere.  The  code  has  yet  to  be  modified  for  dielectric 
surfaces. 

1.2  Problem  Statement 

The  objectives  of  this  research  are  to 

•  modify  the  Wright  Lab  FVTD  code  to  analyze  the  electromagnetic  scattering  from 
complicated  perfect  electric  conductor  (PEC)  three-dimensional  (3-D)  objects 

•  validate  the  characteristic-based  FVTD  formulation  by  using  the  modified  code  to 
analyze  the  scattering  from  the  3-D  objects,  ogive  and  cone-sphere,  and  compare  the 
FVTD  RCS  results  to  MoM  results  and  empirical  data  published  by  the  EMCC. 

Prior  to  the  thesis  research,  the  Wright  Lab  FVTD  code  required  modification  to  analyze 
the  electromagnetic  scattering  from  complicated  three-dimensional  objects  such  as  an  ogive  and 
cone-sphere  or  other  closed-surface,  single-zone  objects.  The  original  FVTD  code  was  limited  to 
calculating  the  RCS  from  a  sphere  and  required  modification  for  closed-surface  three- 
dimensional  objects.  A  modification  changing  the  spherical  coordinate  system  dependency  to  a 
curvilinear  grid  system  which  is  applicable  for  closed-surface  objects  permits  the  code  to  be 
applied  to  complicated  3-D  geometries.  Code  modifications  also  include  the  ability  to  add  grid 
points  near  tips  of  objects  to  ensure  accuracy.  The  incident  wave  specification  in  the  original 
FVTD  code  also  required  modification.  The  original  code  limited  the  specification  of  the 
incident  wave  to  a  sinusoid  with  one  polarization  and  one  direction  of  propagation. 
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Modifications  permit  an  incident  wave  with  options  for  specifying  the  direction,  polarization, 
and  type,  either  a  sinusoid  or  broadband  Gaussian  pulse.  The  incident  wave  modifications  permit 
the  calculation  of  monostatic  and  bistatic  RCS  results  for  different  polarizations.  The  Gaussian 
pulse  provides  scattering  results  for  multiple  frequencies  with  one  simulation.  A  bistatic-to- 
monostatic  approximation  reduces  computer  simulation  times  for  monostatic  RCS  calculations. 

An  analysis  of  the  effectiveness  of  using  FVTD  to  compute  the  electromagnetic 
scattering  from  complicated  objects  is  a  requirement  before  the  Wright  Lab  FVTD  code  and 
algorithm  can  be  considered  a  feasible  computational  design  tool.  The  comparison  of  FVTD 
RCS  calculations  to  a  historically  proven  CEM  method  such  as  the  Method  of  Moments  and  to 
empirical  RCS  data  validates  its  capability  and  accuracy.  This  validation  is  completed  for  the 
EMCC-defined  ogive  and  cone-sphere  in  this  thesis  research. 

1.3  Summary  of  Current  Knowledge 

The  FVTD  formulation  is  characteristic-based.  Chapter  2  gives  an  overview  of  the 
characteristic-based  schemes  used  by  different  researchers.  Dr.  Shang’s  specific  FVTD 
formulation  and  characteristic-based  numerical  algorithm  are  described  in  detail  in  Appendix  A. 
Reportedly,  the  advantages  of  characteristic-based  computational  algorithms  are  [30] 

•  a  naturally  enforced  well-posedness  condition  for  specifying  initial  values 

•  a  windward  spatial  discretization  based  on  the  direction  of  wave  propagation 

•  a  radiation  boundary  condition  (BC)  based  on  a  compatibility  condition 
Specifically,  the  windward  spatial  discretization  imitates  the  physics  of  the  electromagnetic 
scattering.  For  instance,  it  imitates  the  direction  of  wave  propagation.  The  sign  of  the 
eigenvalue  corresponds  to  the  direction  of  propagation.  Forward  differencing  is  used  with  the 
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negative  eigenvalues  and  backward  differencing  with  the  positive  eigenvalues.  The  compatibility 
condition  minimizes  erroneous  numerical  reflections  from  the  outer  boundary  of  the  grid. 

The  process  of  applying  the  characteristic-based  FVTD  formulation  to  electromagnetics 
problems  is  sequential  and  follows  a  logical  procedure.  First,  the  physical  space  surrounding  the 
object  of  interest,  such  as  an  aircraft,  is  discretized  into  finite-volume  cells.  The  space  grid  refers 
to  the  computational  space  containing  the  cells  [51].  The  frequency  of  interest  and  the  electrical 
length  of  the  object  determines  the  number  of  cells  required  in  the  space  grid  [47].  The  finite 
truncated  space  of  a  computational  domain  will  generate  numerical  reflections  in  the  space  and 
produce  erroneous  computations  [37].  Radiation  boundary  conditions,  such  as  the  compatibility 
condition  which  sets  the  incoming  flux  component  equal  to  zero,  are  implemented  at  the  edges  of 
the  truncated  space  and  reduce  the  numerical  reflections  in  the  grid  space. 

Second,  for  use  in  the  FVTD  algorithm,  the  time-domain  Maxwell  partial  differential 
equations  are  placed  in  conservation  form  (See  Chapter  2  or  Appendix  A)  [47].  Maxwell’s 
equations  comprise  a  hyperbolic  system  of  partial  differential  equations  which  can  be  solved 
using  characteristic-based  analysis.  The  Maxwell  equations  in  conservation  form  are  applied  to 
every  finite-volume  cell  in  the  grid  and  solved  numerically  in  the  spatial  and  temporal  domains 
using  one  of  several  techniques.  The  spatial  techniques,  or  flux  evaluation  methods,  are  broadly 
categorized  as  implicit  or  explicit  (See  Chapter  2)  [37].  For  the  time  integration.  Maxwell’s 
equations  in  FVTD  form  are  integrated  using  one  of  several  techniques. 

For  the  implementation  of  the  FVTD  technique,  the  scattered-field  formulation  is  used  as 
opposed  to  the  total-field  formulation.  Only  the  scattered  electric  and  magnetic  fields  and  not  the 
total  fields  propagate  through  the  space  grid.  The  scattered-field  formulation  is  implemented  for 
a  PEC  object  by  setting  the  tangential  electric  field  on  the  scatterer  surface  equal  to  the  negative 
of  the  tangential  incident  field.  The  incident  field  never  appears  in  the  FVTD  calculations. 
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Finally,  the  far-field  RCS  is  computed.  The  scattered  fields  computed  by  FVTD  give 
calculations  in  the  near-field  and  in  the  time-domain.  A  Fourier  transform  permits  frequency 
data  to  be  obtained  from  the  time-domain  data.  A  Green’s-function-based  transformation  allows 
the  far-field  RCS  to  be  easily  calculated  from  the  near-field  frequency  data  [51]. 

1.4  Assumptions 

To  modify  the  Wright  Lab  FVTD  code  for  complicated  objects,  several  assumptions 
were  made.  First,  the  radiation  boundary  condition,  or  compatibility  condition,  implemented  by 
Shang  was  assumed  to  be  sufficient.  The  compatibility  condition  sets  the  incoming  flux 
component  equal  to  zero.  The  boundary  condition  is  only  accurate  if  the  scattered 
electromagnetic  wave  is  parallel  to  one  of  the  transformed  coordinates  [30,  38].  Numerical 
errors  introduced  because  the  propagating  wave  is  not  parallel  to  a  coordinate  axis  were  not 
addressed.  The  first-order  surface  boundary  condition  for  a  PEC  surface  was  also  assumed  to  be 
sufficient  for  the  purposes  of  this  thesis  research. 

Second,  the  mathematical  algorithm  used  by  Shang,  van  Leer’s  kappa  scheme,  was 
assumed  to  be  the  most  adequate  method  for  the  flux  evaluation  in  Maxwell’s  equations.  Shang 
has  studied  numerous  FVTD  numerical  techniques  for  solving  Maxwell’s  equations  and  the 
current  flux  evaluation  technique  is  assumed  to  be  the  most  efficient  and  accurate  [29-41]. 

1.5  Scope 

The  focus  of  this  thesis  research  was  to  modify  the  Wright  Lab  FVTD  code  to  compute 
the  scattering  from  the  ogive  and  cone-sphere  and  validate  the  FVTD  code  by  comparing  FVTD 
RCS  results  to  MoM  and  empirical  RCS  data. 
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Perfect  electric  conductor  (PEC)  surfaces  are  used  for  the  complicated  three-dimensional 
objects.  Dielectric  surfaces  were  not  studied;  therefore,  multi-zoning  which  requires  generating 
different  grids  for  each  dielectric  surface  or  layer,  is  not  implemented  in  the  FVTD  code. 

FVTD  codes  potentially  require  massive  amounts  of  computational  time  and  memory. 
Modifying  the  code  for  speed  optimization  or  parallel  computing  was  not  emphasized  because 
other  researchers  such  as  Blake  are  optimizing  the  code  for  parallel  computing  [6,  7]. 

1.6  Thesis  Organization 

With  Chapter  1  serving  as  the  basic  foundation  for  FVTD  and  the  research  which  needs 
to  be  completed  in  the  area  of  CEM,  Chapter  2  discusses  the  FVTD  formulation  and  the 
numerical  techniques  used  by  various  researchers  such  as  Shang  and  Shankar  to  solve 
electromagnetic  scattering  problems.  Chapter  3  discusses  the  methodology  used  in  completing 
the  research  and  the  code  modifications  that  were  required  to  thoroughly  analyze  the  scattering 
from  the  ogive  and  cone-sphere  using  FVTD.  The  FVTD  RCS  results  for  the  ogive  and  cone- 
sphere  using  the  modified  Wright  Lab  FVTD  code  are  presented  in  Chapter  4  along  with  a 
discussion  of  the  results.  Also  in  Chapter  4,  FVTD  results  are  compared  to  MoM  results  and 
empirical  RCS  data.  Chapter  5  contains  conclusions  on  the  FVTD  research  and  includes 
proposed  areas  of  future  FVTD  research.  The  two  appendices  contain  supplemental  data  on  the 
FVTD  formulation  and  code  development.  In  Appendix  A,  the  FVTD  formulation  and  numerical 
algorithm  implemented  by  Dr.  Shang  in  the  Wright  Lab  FVTD  code  are  discussed.  Appendix  B 
contains  code  listings  and  descriptions  of  the  FVTD  code  modifications. 
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2  Literature  Review 


2.1  Overview 

The  FVTD  computational  technique  is  capable  of  concurrently  solving  the  Euler 
equations  of  fluid  dynamics  and  the  Maxwell  equations  of  electromagnetics.  CFD  has  used  the 
FVTD  technique  since  the  early  1980’s  [48]  to  analyze  the  airflow  about  an  aircraft  or  airfoil  and 
the  technique  has  recently  been  applied  to  CEM.  Several  engineers,  Blake,  Shang,  Shankar, 
Bishop,  Huh,  and  Noack,  are  exploring  and  advancing  the  application  of  the  FVTD  technique  to 
the  Maxwell  equations  of  electromagnetics  [5,  6-7,  15,  22,  29-41,  42-47], 

The  FVTD  technique  follows  a  logical  procedure  from  grid  generation  and  formulating 
the  Maxwell  equations  in  FVTD  form  to  evaluating  the  fields  or  fluxes  through  each  cell  of  the 
grid.  The  computed  scattered-field  results  are  then  transformed  from  the  near-field  to  the  far- 
field.  From  the  far-field  calculations,  RCS  results  are  obtained.  The  following  sections  discuss 
the  FVTD  procedures  for  applying  FVTD  to  electromagnetic  scattering  problems  and  the 
characteristic-based  FVTD  formulations  and  numerical  algorithms  implemented  by  Shang  and 
Shankar.  Appendix  A  describes  in  detail  one  specific  formulation  and  numerical  algorithm  used 
by  Shang  and  implemented  in  the  FORTRAN  code  in  this  thesis  research. 

2.2  Grid  Generation  of  Finite- Volume  Cells 

Computer  simulations  require  that  3-D  geometries  in  a  physical  space  be  accurately 
represented  in  a  computational  domain  [47].  To  use  FVTD,  the  physical  space  surrounding  an 
object  of  interest,  such  as  the  ogive  and  cone-sphere  used  in  this  thesis,  must  be  discretized  into 
volumetric  cells.  The  space  containing  the  finite- volume  cells  is  referred  to  as  the  space  grid 
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[51].  The  frequency  of  interest  and  the  electrical  length  of  the  object  determines  the  number  of 
cells  in  the  grid  [47].  The  grid  can  take  on  several  forms,  either  structured  or  unstructured. 

A  structured  grid  is  defined  by  clearly  distinguishable  coordinate  directions  [7].  Simple 
shapes  or  surfaces  which  can  align  with  the  axes  of  the  three-dimensional  coordinate  system  use 
the  structured  grid.  In  contrast  to  the  structured  grid,  an  unstructured  grid  contains  undefinable 
coordinate  directions  [7].  An  advantage  of  the  unstructured  grid  is  its  ability  to  conform  the  grid 
to  the  surfaces  of  irregular  objects. 

For  characteristic-based  FVTD  formulations,  a  structured  grid  using  curvilinear 
coordinates  is  used  so  the  wave  propagation  is  aligned  closely  with  one  of  the  coordinate  axes 
[34].  The  compatibility  condition  used  for  the  radiation  boundary  condition  is  exact  if  the  wave 
propagation  parallels  a  coordinate  axis.  In  addition,  the  curvilinear  coordinates  permit  higher 
accuracy  in  the  computation  of  the  electric  and  magnetic  scattered  fields.  The  fields  at  the 
centers  of  the  cells  and  the  fluxes  at  the  faces  of  the  cell  are  calculated  using  the  metrics  of  the 
cell  and  the  curvilinear  coordinates  better  represent  the  geometric  features  of  the  object  resulting 
in  higher  accuracy  of  the  scattered  field  computations. 

2.3  Maxwell’s  Equations 

The  time-domain  Maxwell  equations,  in  differential  form,  are  shown  below  and  will  be 
used  in  the  development  of  the  electromagnetic  FVTD  equations: 


Faraday’s  Law: 

VxE  =  -  — 
dt 

(1) 

Ampere’s  Law: 

VxH=^  +  J 
dt 

(2) 

Gauss’s  Electric  Law: 

V-Z)  =  p 

(3) 

Conservation  of  Magnetic  Charge: 

V- jB  =  0 

(4) 

10 


where  E\  Electric  field  strength  vector  (V/m) 

D\  Electric  flux  density  vector  (C/m^) 

H\  Magnetic  field  strength  vector  (A/m) 

B\  Magnetic  flux  density  vector  (Wb/m^  or  T) 

J:  Electric  current  density  vector  (A/m^) 

p:  Electric  charge  density  {Cirri) 

The  constitutive  parameters  relate  the  field  strength  vectors  and  the  flux  density  vectors. 
The  constitutive  parameters,  the  electric  permittivity  and  the  magnetic  permeability,  are  normally 
expressed  as  tensors.  However,  if  the  material  is  linear  and  isotropic,  the  constitutive  parameters 
are  scalars  and  the  constitutive  relations  become 

D  =  eE  and  B  =  pH  (5) 

where  e:  Electric  permittivity  (F/m) 

p;  Magnetic  permeability  (H/m) 

The  four  Maxwell  equations  are  interdependent.  The  two  divergence  equations  can  be 
derived  from  the  two  curl  equations  using  the  conservation  of  charge  relationship 
V  -  J  =  -  3p  /  3t  assuming  the  material  is  linear  and  isotropic.  The  FVTD  calculations  do  not 
use  the  two  divergence  equations,  but  the  equations  can  be  used  as  a  check  on  the  predicted  field 
response  [18]. 

2.4  Maxwell’s  Equations  in  Conservation  Form 

For  use  in  FVTD,  the  two  curl  Maxwell  equations  are  cast  in  conservation  form  [37]. 

The  solution  of  Maxwell’ s  equations  do  not  require  the  conservation  form;  however,  the  form  is 
required  by  the  Euler  equations  to  conserve  physical  properties  such  as  energy,  mass,  and 
momentum  [8].  The  Maxwell  equations  are  also  cast  in  conservation  form  to  take  advantage  of 
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the  same  computational  technique  used  to  solve  the  Euler  equations.  To  place  the  two  curl 
Maxwell  equations  in  conservation  form,  the  curl  operations  are  carried  out  and  the  constitutive 
parameters  are  implemented.  The  result  is  given  by  [37] 

du  dF  8G  an  , 


where 


Equation  (6)  is  a  system  of  six  linear  equations.  U  is  the  independent  variable  and  the  F,  G,  and 
H  flux  vectors  are  the  dependent  variables.  The  equations  are  not  linearly  independent; 
therefore,  a  characteristic-based  technique  is  needed  to  uncouple  the  six  equations. 


2.5  Coordinate  Transformation 

To  analyze  the  scattering  of  various  objects,  including  the  ogive  and  cone-sphere 
analyzed  in  this  thesis,  a  curvilinear  coordinate  transformation  is  required  [37].  A  curvilinear 
structured  grid  minimizes  the  errors  introduced  in  the  cell  metrics  and  the  flux  calculations.  The 
order  of  accuracy  will  be  below  the  formal  order  of  accuracy  if  a  poor  curvilinear  grid  is 
generated  [34].  The  coordinate  transformation  converts  the  Cartesian  coordinates  to  curvilinear 
coordinates.  The  transformation  defines  a  one-to-one  relationship  between  two  sets  of  temporal 
and  spatially  dependent  variables.  The  variables  used  are  r],  and  ^  and  are  all  functions  of  x,  y, 
and  z.  Equation  (6),  after  a  coordinate  transformation,  becomes  [37,  38] 
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(7) 


H  = 


V 


V  is  the  Jacobian  of  the  coordinate  transformation  and  is  equal  to  [34] 
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A  coordinate  transformation  closely  aligns  the  direction  of  propagation  of  the  scattered 
electric  and  magnetic  fields  to  the  coordinate  directions  [37].  The  radiation  boundary  condition 
produces  fewer  numerical  reflections  if  the  coordinate  axis  closely  parallels  the  direction  of 
propagation  of  the  scattered  fields. 


2.6  Finite-Volume  Formulation 

Equation  (7)  is  applied  to  every  finite-volume  cell  in  the  grid.  An  integration  is 
performed  over  each  finite- volume  cell: 
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(9) 


fff  dU  fff  dF  dG  dff 

JJJv  at  ^  ari  ac 


dV 


The  divergence  theorem  is  then  applied  to  the  second  integral: 


IJJvf  =-1/1/'“^  ™ 

where  n:  Unit  vector  normal  to  the  surface  (4,  Ti,  and  4  for  F,  G,  and  H,  respectively) 

S:  Closed  surface  bounding  the  finite  volume  (m^) 

Equation  (10)  is  the  expression  for  a  generic  FVTD  formulation.  The  unknown 
components  of  the  U  vector  are  the  magnetic  and  electric  flux  densities.  The  vectors 
F,  G,  and  H  are  the  flux  vectors  and  can  be  expressed  in  terms  of  the  magnetic  and  flux 
densities.  A  multitude  of  techniques  are  used  to  solve  Equation  (10)  and  differentiate  the  myriad 
of  FVTD  numerical  algorithms. 


2.7  Boundary  Conditions 

Realistically,  radar  waves  scattered  from  an  object  travel  away  from  the  object  and  are 
not  reflected  back  to  the  target.  However,  the  truncated  space  of  a  computational  domain  will 
generate  numerical  reflections  in  the  space  and  produce  erroneous  computations  [45]. 
Implementing  a  radiation  boundary  condition  (BC)  reduces  the  numerical  reflections  in  the  grid 
space. 

Shankar  and  Shang  use  a  first-order  accurate  radiation  BC  [30,  45],  For  precise 
calculations,  first-order  accuracy  is  not  acceptable  and  higher-order  BCs  would  ensure  higher 
accuracy.  Shang  and  Shankar  use  a  compatibility  condition  in  which  the  incoming  flux 
component  is  set  to  zero  at  the  boundary  [37].  For  the  compatibility  condition,  the  fields 
traveling  perpendicular  to  the  boundary  are  not  reflected.  For  example,  in  the  case  of  the 
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propagation  of  a  wave  from  a  dipole,  the  BC  is  exact  since  the  wave  travels  along  the  radial 
coordinate  direction.  However,  numerical  errors  can  result  if  the  wave  is  not  traveling 
perpendicular  to  the  boundary.  The  coordinate  transformation  discussed  previously  increases  the 
component  of  the  wave  traveling  perpendicular  to  the  outer  boundary  [37]. 

A  surface  boundary  condition  is  implemented  on  the  surface  of  PEC  scatterers.  The  BC 
sets  the  tangential  electric  field  equal  to  zero  and  the  normal  component  of  the  magnetic  flux 
density  equal  to  zero  [7].  Details  for  the  PEC  BCs  implemented  for  this  thesis  research  can  be 
found  in  Appendix  A.  Boundary  conditions  are  also  required  for  dielectric  interfaces  and  details 
for  these  BCs  can  be  found  in  the  Shankar  references  [42-47]. 

2.8  Flux  Evaluation  and  Time  Integration 

The  flux  vectors  in  Equation  (10)  can  be  evaluated  numerically  using  one  of  several 
techniques  which  can  be  broadly  categorized  as  explicit  or  implicit  [30,  47].  Explicit  numerical 
expressions  place  the  independent  and  dependent  variables  on  different  sides  of  the  equation. 
Implicit  expressions  are  recognized  by  the  intermixing  of  the  dependent  and  independent 
variables  on  each  side  of  the  equation. 

The  techniques  implemented  by  Shang  and  Shankar  are  explicit  characteristic-based 
schemes.  The  objective  of  the  characteristic-based  numerical  procedures  is  to  achieve  the 
Riemann  approximation  to  a  three-dimensional  problem  in  each  spatial  dimension  [30].  The 
three-dimensional  Maxwell  equations  can  then  be  solved  in  each  dimension  sequentially. 

Shang  has  studied  and  applied  several  characteristic-based  implicit  and  explicit 
techniques.  He  is  presently  using  an  explicit  van  Leer’s  kappa  scheme  in  which  the  flux  on  a 
surface  of  a  cell  is  extrapolated  from  data  of  adjacent  cell  centers  [34].  The  scheme  is  referred 
to  as  a  Monotone  Upstream-Centered  Scheme  for  Conservation  Laws  (MUSCL)  and  is  a 
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windward  approach  that  considers  the  direction  of  wave  propagation.  The  MUSCL  approach 
produces  various  orders  of  accuracy.  The  approach  used  in  this  thesis  research  produces  results 
which  are  third-order  accurate.  The  details  of  van  Leer’s  kappa  scheme  are  discussed  in 
Appendix  A. 

A  flux- vector  splitting  algorithm  developed  by  Steger  and  Warming  [48]  is  used  to 
calculate  the  fluxes  from  the  independent  variable  U .  The  incoming  and  outgoing 
electromagnetic  waves  are  split  based  on  the  positive  and  negative  sign  of  the  characteristic, 
hence,  the  name  split-flux  vectors. 

Shankar  uses  an  explicit  Lax-Wendroff  upwind  scheme  that  is  characteristic-based.  The 
scheme  uses  a  predictor  and  a  corrector  step.  The  predictor  step  results  in  a  first-order  accurate 
solution.  The  corrector  step  increases  the  accuracy  of  the  solution  to  second-order  [42]. 

Equation  (10),  in  the  temporal  or  time-stepping  domain,  can  be  solved  using  several 
techniques,  just  as  in  the  spatial  domain.  Shankar  uses  the  same  second-order  Lax-Wendroff 
upwind  scheme.  Shang  uses  a  Runge-Kutta  family  of  single-step  multi-stage  procedures  [32,  38] 
which  gives  varying  degrees  of  accuracy.  For  example,  with  van  Leer’s  kappa  scheme  for  the 
flux  evaluation,  he  uses  either  a  two-stage  or  four-stage  Runge-Kutta  method  that  produces 
second-order  accuracy  and  fourth-order  accuracy,  respectively  [34]. 

The  higher-order  accurate  Wright  Lab  code  can  potentially  reduce  the  required  number 
of  grid  points  per  wavelength  permitting  the  computation  of  the  electromagnetic  scattering  from 
electrically  larger  objects  and  reducing  the  computer  simulation  time  [36, 40]. 

2.9  Green’s-Function-Based  Near-to-Far  Field  Transformation 

The  spatial  and  time  integration  of  Equation  (10)  gives  results  in  the  near-field.  Various 
calculations,  such  as  the  RCS,  are  far-field  results.  The  FVTD  grid  is  in  the  near  field;  however. 
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Green’  s-function-based  transformations  allow  the  scattered  fields  in  the  far-field  to  be  easily 
calculated  from  the  near-field  results  subsequent  to  a  Fourier  transform  [51],  The  near-to-far 
field  transformation  permits  the  FVTD  grid  to  be  truncated  in  the  near-field  and  does  not  have  to 
extend  out  to  the  far-field  to  obtain  the  RCS  data. 

The  far-field  results  are  obtained  by  creating  a  virtual  surface  around  the  object.  This 
surface  does  not  have  to  conform  to  the  object.  An  imaginary  surface  in  the  FVTD  grid  space 
can  serve  as  a  virtual  surface.  The  surface  equivalence  theorem  is  applied  to  the  surface  to 
obtain  the  equivalent  time-harmonic  electric  and  magnetic  currents  and  charges.  The  currents 
and  charges  on  the  virtual  surface  are  then  weighted  by  a  free-space  Green’s  function  to  obtain 
the  far-field  E  and  H  fields  [51].  The  far-field  results,  such  as  the  RCS,  are  easily  calculated 
from  the  far-field  scattered  E  and  H  fields. 

2.10  Applications 

As  mentioned  previously,  FVTD  is  relatively  new  to  CEM.  Shankar,  Shang,  Blake,  and 
others  have  applied  the  FVTD  codes  to  analyzing  the  electromagnetic  scattering  from  various 
objects.  Shankar  has  analyzed  a  two-dimensional  PEC  and  dielectric  circular  cylinder  [44].  Also 
in  2-D,  Shankar  has  analyzed  a  simple  engine  inlet,  ogive,  and  airfoil  [44].  In  addition,  Shankar 
has  completed  extensive  3-D  research  in  which  he  has  analyzed  a  PEC  sphere  and  missile.  He 
has  also  worked  with  frequency-dependent  materials  [44].  Other  applications  can  be  found  in 
Shankar’s  contract  report  [42]. 

Shang  has  used  his  various  FVTD  codes  to  analyze  the  scattering  from  a  sphere, 
propagation  from  a  dipole,  and  the  propagation  of  a  wave  through  a  waveguide.  Blake  has 
rewritten  Shang’ s  FVTD  vectorized  code  for  parallel  computing  machines  and  has  used  his 
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modified  code  to  study  the  propagation  of  an  electromagnetic  wave  in  a  waveguide  and  the 
scattering  from  a  single  sphere  and  dual  spheres  [6,  7]. 

2.11  Summary 

Computational  techniques  which  permit  simultaneous  trade-offs  between  the 
electromagnetic  and  aerodynamic  disciplines  would  greatly  improve  the  efficiency  of  the 
engineering  design  process  for  low-observable  aircraft.  The  FVTD  computational  technique  is 
capable  of  solving  the  Euler  equations  of  fluid  dynamics  and  the  Maxwell  equations  of 
electromagnetics.  The  FVTD  technique,  historically  proven  for  fluid  dynamics,  has  excellent 
potential  in  computational  electromagnetics.  Several  engineers,  such  as  Shang,  Shankar,  and 
Blake,  are  advancing  the  application  of  the  FVTD  technique  to  solving  the  Maxwell  equations  of 
electromagnetics . 

The  FVTD  technique  follows  a  logical  process  for  calculating  electromagnetic  data  for  a 
scattering  object.  A  curvilinear  grid  is  generated  about  an  object  to  ensure  accurate  results  from 
the  FVTD  formulation.  The  Maxwell  equations  are  placed  in  conservation  form  and  then  solved 
in  the  spatial  and  temporal  domains  using  a  technique  dependent  on  the  desired  accuracy. 
Numerical  reflections  in  the  domain  are  reduced  by  implementing  outer  boundary  conditions. 
Depending  on  the  desired  results,  transformations  can  be  applied  to  the  computational  results. 
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3  Methodology 


3.1  Overview 

For  the  thesis  research,  a  disciplined  approach  was  taken  to  modify  the  code  and  validate 
the  FVTD  formulation  by  studying  the  electromagnetic  scattering  from  an  ogive  and  cone-sphere. 
Figure  1  illustrates  the  basic  approach  used  to  complete  the  FVTD  research.  By  performing 
several  tasks  concurrently,  more  research  was  completed  successfully  and  efficiently.  The 
columns  in  the  block  diagram  correspond  to  the  tasks  which  were  performed  in  parallel  and 
include  the  broad  areas  of  FVTD  code,  comparisons/benchmarks,  and  computer  support. 


Figure  1:  FVTD  Thesis  Research  Block  Diagram 


The  first  task,  FVTD  code,  was  the  focus  of  the  thesis  research.  The  FVTD  code  was 
first  used  to  solve  for  the  scattering  of  a  PEC  sphere  to  become  familiar  with  the  software,  the 
FVTD  formulation,  grid  generation,  and  boundary  conditions.  Analytical  solutions  for  the 
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scattering  of  a  sphere  were  used  as  benchmarks  and  for  verifying  code  modifications  to  ensure 
accurate  changes.  Shang  has  published  results  for  the  sphere  using  the  original  FVTD  code  [38], 
The  reference  includes  comparisons  to  analytical  solutions. 

After  becoming  familiar  with  the  code,  it  was  modified  to  obtain  electromagnetic 
scattering  results  from  the  PEC  surfaces  of  an  ogive  and  cone-sphere.  The  modified  code  can  be 
used  for  other  closed-surface  perfect  electric  conductor  (PEC)  3-D  objects  if  an  appropriate  grid 
is  generated.  Appendix  B  contains  the  detailed  requirements  for  the  grid.  The  Rockwell 
International  Science  Center  generated  the  grid  for  the  ogive  and  the  Flight  Dynamics  Directorate 
of  Wright  Labs  generated  the  grid  for  the  cone-sphere.  In  addition  to  modifying  the  code  for  a 
generic  3-D  object,  the  option  of  specifying  the  direction,  polarization,  and  type  of  the  incident 
wave  was  programmed.  Also  programmed  was  an  RCS  convergence  check  used  with  a  sinusoid 
incident  wave  which  ends  the  simulation  after  the  transients  diminish.  A  threshold  check,  used 
with  a  Gaussian  pulse  incident  wave,  ends  the  simulation  once  the  scattered  field  is  below  a  pre¬ 
determined  threshold.  A  bistatic-to-monostatic  approximation  saves  computer  simulation  time 
for  monostatic  computations. 

To  complete  the  research,  another  CEM  technique  and  experimental  measurements  were 
used  as  benchmarks.  Column  two  in  Figure  1  corresponds  to  the  code  and  experimental  results 
used  as  comparisons.  The  Method  of  Moments,  considered  as  a  very  accurate  CEM  code,  is  used 
as  a  benchmark.  The  body-of-revolution  Moment  Method  code,  CICERO  [55],  was  used  to 
obtain  bistatic  and  monostatic  RCS  results  for  the  ogive  and  cone-sphere.  Experimental 
measurements  published  by  the  Electromagnetic  Code  Consortium  (EMCC)  are  benchmarks  for 
the  monostatic  calculations  [55].  The  empirical  results  are  excellent  data  to  use  as  benchmarks 
and  validate  developmental  CEM  codes  such  as  the  FVTD  code. 
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The  third  task,  computer  support,  was  critical  in  the  completion  of  the  research.  The  Air 
Force  Institute  of  Technology  (AFIT)  machines,  the  SUN  Sparc  20  and  DEC  Alpha  machines, 
were  used  to  assess  code  modifications.  A  Cray  90  at  the  USAE  Waterways  Experiment  Station 
(CEWES)  high  performance  computing  center  (HPCC)  in  Vicksburg,  Mississippi  was  used  to 
complete  the  majority  of  the  computer  simulations  for  the  ogive.  The  SP-2  at  the  Maui  HPCC  in 
Maui,  Hawaii,  was  used  to  complete  the  monostatic  tests  for  the  cone-sphere.  The  Cray  90  and 
SP-2  were  invaluable  in  completing  the  research. 

3.2  FVTD  FORTRAN  Code  Modifications 

The  original  FVTD  FORTRAN  code  calculated  the  electromagnetic  scattering  from  a 
sphere.  The  code  modifications  permit  the  calculation  of  the  scattering  from  other  closed- 
surface,  single-zone,  3-D  objects.  The  capability  was  added  to  read  grid  files  generated  by 
computer-aided  design  (CAD)  software  programs.  Options  for  changing  the  size  of  the  grid, 
such  as  the  addition  of  grid  points  at  diffraction  points,  were  implemented  to  improve  accuracy 
of  the  scattered-field  computations.  The  original  code  specified  the  incident  wave  to  propagate 
from  only  one  direction  with  one  polarization.  The  direction  and  type  of  incident  wave  was 
added  along  with  an  RCS  convergence  check  and  a  threshold  check  for  the  scattered  field.  A 
bistatic-to-monostatic  approximation  obtains  monostatic  results  from  bistatic  results.  Appendix 
B  contains  the  code  and  subroutines  written  for  the  code  modifications  and  related  functions. 

3.2.1  Geometry:  Ogive 

The  ogive  is  a  very  common  test  body  for  code  validation  and  is  a  classic  low-observable 
shaped  body.  The  ogive  is  a  body-of-revolution  formed  by  rotating  a  convex  arc  around  a  chord 
[27].  No  analytical  solution  for  the  electromagnetic  scattering  from  an  ogive  is  available; 
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therefore,  MoM  results  and  experimental  data  are  used  as  benchmarks.  The  EMCC-defined 
ogive  is  10  inches  (0.254  m)  long,  has  a  maximum  radius  of  1  inch  (0.0254  m),  and  a  half-angle 
of  22.62  degrees  [55].  The  single  ogive  with  a  metallic  surface  is  described  mathematically  as 


f(z)  =  -cos(22.62°)-l- 

_  f(z)cos\|/ 
l-cos(22.62°) 

_  f(z)sin\j/ 
l-cos(22.62°) 

for  -5.0  <  z  <  5.0  inches 
-71  <  'P  <  71  radians 


(11) 

(12) 

(13) 


y  (m) 


Figure  2:  Mesh  Plot  of  EMCC-Defined  Ogive 
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Figure  2  is  a  mesh  plot  of  the  single  ogive.  Note  that  the  ogive  is  approximately  one  wavelength 
long  at  1.18  GHz  [55].  The  units  in  the  plot  are  in  meters. 


3.2.2  Geometry:  Cone-sphere 

The  cone-sphere  is  another  common  RCS  test  body  which  is  also  a  body-of-revolution. 
The  EMCC-defined  cone-sphere  has  a  half-angle  of  7.0  degrees,  radius  of  2.947  inches  (0.07485 


y(m) 


Figure  3:  Mesh  Plot  of  EMCC-Defined  Cone-Sphere 


m),  and  length  of  27.127  inches  (0.6890  m)  [55].  The  metallic  surface  of  the  cone  portion  of  the 
cone-sphere  is  described  mathematically  as 
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10.10 


X  =  0.12279cos(\|/)(z  + 23.821) 


(14) 


y  =  0.12279sin(\{f)(z  +  23.821) 

for  -23.821  <  z  <  0.0  inches 
-TZ<^  <%  radians 

The  surface  of  the  sphere  portion  of  the  cone-sphere  is  described  mathematically  as 


X  =  2.947  cos(\|/ ). 


z- 0.359  Y 
2.947  > 


y  =  2.947  sin(v|r). 


z- 0.359  Y 
2.947  J 


(15) 


(16) 


(17) 


for  0.0  <  z  <  3.306  inches 
-7t  <  'P  <  71  radians 

Figure  3  is  a  mesh  plot  of  the  cone-sphere.  The  cone-sphere  is  approximately  two  wavelengths 
long  at  0.869  GHz  [55].  The  units  in  the  plot  are  in  meters. 


3.2.3  Grid  Files 

A  grid  file  for  the  ogive  was  obtained  from  NASA  but  was  originally  generated  by  the 
Rockwell  International  Science  Center.  Slices  of  the  ogive  grid  are  shown  in  Figures  4  and  5. 
The  dimensions  of  the  ogive  equal  the  size  of  the  EMCC-defined  ogive.  Figure  4  is  a  slice  of  the 
grid  in  the  yz  plane.  The  original  grid  size  is  (10-121-30)  in  spherical  coordinates  (R,0,(t)).  As 
seen,  the  radial  lines  of  the  grid  are  approximately  perpendicular  to  the  surface  of  the  ogive.  This 
characteristic  will  produce  more  accurate  results  from  the  first-order  surface  boundary  condition 
(See  Appendix  A).  The  spacing  of  the  cells  in  the  radial  direction  also  increases  with  an  increase 
in  R.  The  larger  grid  spacing  with  an  increases  in  R  minimizes  the  computer  simulation  time  due 
to  the  fewer  number  of  cells  but  retains  the  numerical  accuracy  of  a  finer  grid  near  the  surface. 
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Figure  5:  Slice  of  the  Ogive  Grid  ( 10-121-30)  in  the  xy  Plane 
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Figure  7:  Slice  of  the  Cone-Sphere  Grid  (50-73-45)  in  the  xy  Plane 
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The  original  cone-sphere  grid  was  generated  by  the  Flight  Dynamics  directorate  of 
Wright  Laboratory  using  GRIDGEN,  a  common  CFD  CAD  package.  The  size  of  the  original 
grid  is  (50-73-45)  in  spherical  coordinates  (R,0,(|)).  The  grid  has  the  same  characteristics  as  the 
ogive  grid  such  as  the  increase  in  the  spacing  of  the  cells  in  the  radial  direction  and  tightly 
packed  cells  at  the  surface.  A  hyperbolic  tangent  distribution  was  used  to  generate  the  grid 
spacing  in  the  radial  direction.  The  radial  lines  are  also  approximately  perpendicular  to  the 
surface,  similar  to  the  ogive  grid.  Slices  of  the  cone-sphere  grid  are  shown  in  Figures  6  and  7. 
Figure  6  is  a  slice  of  the  grid  in  the  xz  plane  and  Figure  7  is  a  slice  in  the  xy  plane.  The 
dimensions,  in  meters,  of  the  cone-sphere  match  the  size  of  the  EMCC-defined  cone-sphere. 

3.2.4  Grid  Modifications 

Simple  calculations  included  in  the  code  change  the  size  of  the  original  grid,  grid 
spacing,  and  the  distance  of  the  outer  boundary  from  the  scatterer.  For  the  ogive  and  cone-sphere 
grids,  the  number  of  grid  points  in  the  phi  direction  are  changed  by  finding  the  length  of  the 
radial.  Then  the  desired  spacing  is  calculated  and  (x,  y)  coordinates  are  generated  based  on  the 
same  radial  distance,  desired  spacing,  and  number  of  grid  points. 

For  each  grid,  to  change  the  number  of  grid  points  in  the  radial  direction  and  distance  to 
the  outer  boundary,  the  following  method  is  used: 

R  =  (l-a)R„-i-ccRi  (18) 

where  R:  Distance  to  the  new  coordinate 

Rq:  Distance  to  outer  radial  coordinate  at  the  edge  of  grid 

Rj!  Distance  to  inner  radial  coordinate  at  the  surface  of  the  scatterer 

a:  Scaling  factor  or  increment 
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The  sizes  and  densities  of  the  grids  were  changed  to  analyze  the  effects  of  grid  size  and 
the  number  of  grid  points  per  wavelength.  Grid  points  are  added  by  averaging  the  grid 
coordinates  of  the  surrounding  points.  Grid  points  can  also  be  removed.  A  listing  of  the  code 
containing  the  methods  for  changing  the  grid  sizes  for  the  ogive  and  cone-sphere  is  given  in 
Appendix  B. 

3.2.5  Direction  and  Polarization  of  the  Incident  Wave 

The  incident  wave’s  direction  of  propagation  and  polarization  are  specified  as  shown  in 
Figure  8.  The  direction  of  the  incident  wave,  r,  is  specified  by  the  spherical  coordinates,  6  and  (|). 


The  direction  is  specified  using  the  angle  from  which  the  incident  field  is  propagating.  This 
angle  should  not  be  confused  with  the  angle  of  the  propagation  vector,  k.  The  polarization  is 
specified  with  Ee  and  E,j,  [18].  The  magnitude  of  the  polarization  components  is  unity  and  the 
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amplitude  of  the  wave  is  specified  separately.  The  displacement  unit  vector,  Fd,  is  used  to  define 
the  relative  spatial  delay  for  each  component  of  the  incident  field.  The  vector,  R',  is  used  to 
calculate  the  value  of  the  incident  field  at  each  location  on  the  surface  of  the  scatterer.  Appendix 
B  contains  more  details  on  specifying  the  incident  field. 

3.2.6  Incident  Wave  Type 

3.2.6.1  Sinusoid 

The  sinusoid  incident  wave,  sin(co(t+(R'TD)/c)),  is  specified  easily  at  one  frequency. 

The  sinusoidal  wave  provides  bistatic  RCS  results  for  one  frequency.  At  the  beginning  of  each 
time  step  for  the  scattered-field  formulation,  the  tangential  scattered  electric  field  is  specified  as 
the  negative  of  the  sinusoid  at  the  scatterer  surface.  The  amplitude  of  the  sinusoid  at  one  time 
step  will  depend  on  the  location  of  the  grid  point  on  the  surface.  The  spatial  delay  is  found  by 
taking  the  dot  product  of  R'  and  the  displacement  unit  vector,  Fd.  The  initial  application  of  the 
sinusoidal  wave  at  the  surface  generates  transients  which  must  attenuate  before  accurate 
frequency  data  can  be  obtained.  The  time  for  the  transients  to  die  out  is  discussed  in  Chapter  4. 
The  frequency  data  for  the  near-to-field  transformation  are  taken  after  the  transients  diminish. 


3.2.6.2  Gaussian  Pulse 

An  advantage  of  the  time-domain  analysis  is  the  capability  to  obtain  multiple  frequency 
results  simultaneously.  A  Gaussian  pulse  is  ideal  to  use  as  an  incident  wave  for  multiple 
frequency  analysis.  The  Gaussian  pulse  is  specified  as  [18] 


g(t)  =  exp 


t-t(R^-Fo)/c-tp-t„Y 


T 


) 


(19) 
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where  c:  Velocity  of  the  incident  wave  in  free  space 

g(t):  Gaussian  pulse  with  an  amplitude  of  unity 

to:  Center  of  Gaussian  pulse,  mean 

to:  Delay  of  Gaussian  pulse  due  to  maximum  spatial  delay  in  the  direction  of  the 

incident  wave  (IRdI/c) 

T:  Duration  of  the  Gaussian  pulse 

The  parameters,  T  and  to,  of  the  Gaussian  pulse  can  be  specified  to  obtain  accurate 
results  for  a  desired  bandwidth  of  frequencies  [18].  T  is  the  duration  of  the  pulse.  The 
parameter  to  determines  the  delay  of  the  center  of  the  pulse  and  the  amplitude  of  the  pulse  at 
truncation.  The  pulses  used  for  the  thesis  research  are  tmncated  140  dB  from  the  peak  of  the 
pulse  [18].  The  location  of  the  tmncation  point  is  approximately  5.7  standard  deviations  from 
the  center  of  the  pulse.  Two  Gaussian  pulses  are  illustrated  in  Figure  9.  The  wider  pulse  has  a 
duration  of  0. 1278  nsec,  delay  of  0.5 1 1 1  nsec,  and  a  bandwidth  (BW)  of  10  GHz.  The  narrower 
pulse  has  a  width  of  0.07099  nsec,  delay  of  0.2840  nsec,  and  a  BW  of  18  GHz.  The  frequency 
spectra  of  the  same  pulses  are  shown  in  Figure  10.  The  useful  range  of  frequencies  is  taken  to  be 
roughly  one- third  of  the  bandwidth  shown  [18].  As  seen  in  Figure  10,  sufficient  amplitude  is 
available  for  that  range  of  frequencies.  Appendix  B  contains  a  reference  for  specifying  the 
parameters  for  a  Gaussian  pulse  incident  wave. 

The  use  of  either  a  Gaussian  pulse  or  a  sinusoid  wave  for  the  incident  field  is  dependent 
on  the  results  desired.  The  simulation  for  a  Gaussian  pulse  is  usually  longer  due  to  the  time  that 
it  takes  the  pulse  to  travel  through  the  computational  space  and  for  the  scattered  field  to  diminish. 
Therefore,  the  Gaussian  pulse  is  advantageous  only  if  the  time  required  for  the  simulation  is  less 
than  the  total  time  required  to  complete  tests  for  individual  sinusoid  cases. 
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Figure  9:  Gaussian  Pulses  with  Different  Bandwidths 


Figure  10:  Frequency  Spectrum  of  Gaussian  Pulses  Shown  in  Figure  9 
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3.2.7  Scattered-Field  Checks 


Figure  1 1  shows  the  decision  flowchart  for  the  scattered-field  checks  added  to  the  FVTD 
code.  The  checks  ensure  that  accurate  frequency  data  is  taken  before  the  simulation  ends.  The 
checks  occur  at  the  end  of  a  time  step  or  a  period.  If  the  checks  result  in  a  “Yes,”  the  simulation 
ends  and  the  RCS  values  are  calculated.  If  the  result  is  “No”,  the  time  loop  is  entered  again  to 
calculate  the  fluxes  for  another  time  step.  If  the  checks  are  not  enabled,  the  Fourier  transform 
parameters  included  in  the  input  file  end  the  simulation  (See  Appendix  B). 


Figure  11:  Flowchart  for  RCS  Convergence  and  Threshold  Checks 
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3.2.7.1  RCS  Convergence 

When  initially  specified,  the  sinusoid  incident  wave  generates  transients  that  must  die 
out  before  accurate  results  can  be  obtained.  To  ensure  convergence  for  the  RCS  values  when  a 
sinusoidal  incident  wave  is  implemented,  a  convergence  check  was  programmed  in  the  code. 

The  check  calculates  the  RCS  value  at  ten  different  viewing  angles.  If  the  RCS  values  are  within 
0. 1  dB  of  the  RCS  values  from  the  previous  time  step,  convergence  has  been  reached.  The 
bistatic  RCS  values  are  calculated  from  the  frequency  data  taken  for  one  period  after 
convergence.  Appendix  B  contains  the  code  listing  for  the  RCS  convergence  check. 

3.2.7.2  Threshold  Check 

A  threshold  check  for  the  scattered  field  resulting  from  a  Gaussian  pulse  incident  wave 
was  added  to  the  code.  The  check  ends  the  simulation  if  the  amplitude  of  the  scattered  field  is 
less  than  140  dB  of  the  maximum  amplitude  of  the  incident  wave.  The  check  is  performed  by 
sampling  the  scattered  field  one  cell  away  from  the  scatterer  surface.  The  virtual  surface  which 
is  used  to  calculate  equivalent  currents  and  the  far-field  RCS  results  is  also  located  one  cell  away 
from  the  target  surface.  Sampling  scattered-field  amplitudes  on  the  virtual  surface  ensures  that 
the  frequency  data  is  accurate  for  the  RCS  calculations. 

3.2.8  Bistatic-to-Monostatic  Approximation 

One  simulation  of  the  FVTD  code  gives  bistatic  data  at  721  receiving  locations.  The 
angles  are  spaced  0.25°  degrees  apart,  providing  a  180°  sweep  of  bistatic  RCS  values.  A 
monostatic  calculation  is  obtained  for  only  one  location.  To  obtain  monostatic  results  for  the 
entire  1 80°  sweep,  a  monostatic  FVTD  test  must  be  performed  for  every  viewing  location  which 
would  require  721  computer  simulations.  Obviously,  this  would  require  an  enormous  amount  of 
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computer  time  and  effort.  A  bistatic-to-monostatic  approximation  provides  the  entire  monostatic 
sweep  and  decreases  the  number  of  monostatic  tests  required  [28], 

The  bistatic-to-monostatic  approximation  is  expressed  as 

RCSm(0.5-  a)  =  RCSgCa)  (20) 

where  RCSb:  Calculated  bistatic  RCS  value 

RCSm:  Approximate  monostatic  RCS  value 
a:  Angle  at  which  bistatic  value  is  calculated 

The  approximation  originated  when  researchers  at  RCS  test  ranges  observed  that  experimental 
monostatic  RCS  values  are  equal  to  bistatic  values  at  twice  the  receiving  angle  [28].  The  angle 
must  be  small,  usually  less  than  or  equal  to  5°. 

The  monostatic  calculations  in  Chapter  4  are  obtained  by  performing  monostatic 
simulations  every  10°.  The  bistatic-to-monostatic  approximation  uses  the  bistatic  data  at  angles 
on  each  side  of  the  incident  angle  to  obtain  the  monostatic  data.  For  example,  the  bistatic  data 
from  0.50°  to  10°  is  used  to  calculate  the  monostatic  data  from  0.25°  to  5°.  As  will  be  seen  in  the 
monostatic  plots,  the  approximation  is  sufficient  for  obtaining  reasonable  RCS  results. 

3.3  FVTD  Computational  Issues 

3.3.1  Scattered-Field  Formulation 

The  FVTD  code  uses  a  scattered-field  formulation  in  which  only  the  scattered  field 
propagates  through  the  computational  space.  The  scattered-field  formulation,  as  opposed  to  the 
total-field/scattered  field  formulation,  avoids  the  numerical  dispersion  and  dissipation  of  the 
incident  wave  as  it  propagates  through  the  grid  [52].  Also,  the  scattered  field  is  analytically 
specified  at  the  scatterer  surface.  A  disadvantage  of  the  scattered-field  formulation  is  that  it 
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requires  the  calculation  of  the  scattered  field  at  every  time  step  at  every  point  on  the  surface  of 
the  scatterer. 

The  scattered-field  formulation  is  implemented  for  a  PEC  object  by  specifying  the 
scattered  tangential  electric  field  as  the  negative  of  the  incident  field  at  the  surface.  The  total 
tangential  electric  field  is  zero  at  the  surface  of  the  PEC  object:  E^=E‘+E*  =  0.  Therefore, 
the  tangential  scattered  electric  field  is  equal  to  the  negative  of  the  incident  electric  field  at  the 
surface:  E"=-E‘. 

The  scattered  field  is  specified  at  the  boundary  of  the  scatterer.  Because  of  the  unknovi^n 
electric  and  magnetic  fields  required  at  the  scatterer  surface  for  the  FVTD  formulation,  an 

extrapolated  boundary  condition  is  used  for  the  magnetic  field.  See  Appendix  A  for  the  details 

‘ 

on  the  extrapolated  boundary  condition. 


3.3.2  Stability 

The  stability  of  the  explicit  FVTD  formulation  is  based  on  the  relationship  between  the 
spatial  discretization  and  the  time  step.  The  time  step  is  directly  proportional  to  the  smallest  grid 
spacing  and  is  inversely  proportional  to  the  wave  velocity.  The  Courant-Friedrichs-Lewy  (CEL) 
number  determines  the  stability  of  the  algorithm  [38].  For  example,  for  a  three-dimensional 
problem,  the  CEL  number  has  a  maximum  of  1.74.  The  stability  condition  is  [51] 


At< 


CFL/c _ 

I  r~ 

H - 7  3 - ~ 

(Ay)'  (Az)' 


{21) 


where  CEL: 
At: 


Courant-Friedrichs-Lewy  number 
Time  step 


Ax,  Ay,  Az:  Spatial  increments 


The  stability  of  a  system  is  also  related  to  the  eigenvalues  of  the  system  of  equations.  An 
advantage  of  the  characteristic-based  formulation  is  that  it  addresses  the  fundamental  issues  of 
well-posedness  and  numerical  stability  of  a  system  of  hyperbolic  differential  equations. 

Several  disadvantages  of  the  explicit  FVTD  formulation  are  due  to  the  stability  criteria. 
First,  surface  characteristics  such  as  tips  and  edges  produce  converging  grid  lines  resulting  in 
extremely  small  cells.  The  small  cell  sizes  result  in  a  very  small  time  step  which  increases 
simulation  times.  Also,  for  higher  frequencies,  the  wavelength  decreases  and  the  grid  spacing 
decreases  to  maintain  the  same  number  of  cells/A,  increasing  the  computation  time. 

3.3.3  Numerical  Dispersion 

The  discretized  space  causes  unwanted  dispersion  of  the  scattered  field.  Numerical 
dispersion  varies  with  the  direction  of  propagation,  frequency,  and  variation  of  the  cell  size  [52]. 
The  amount  of  dispersion  increases  as  the  electrical  size  of  the  object  increases.  To  minimize  the 
effects  of  the  numerical  dispersion  on  the  RCS  results,  the  virtual  surface  used  for  calculating 
equivalent  currents  is  located  close  to  the  scatterer,  usually  one  cell  away  from  the  surface. 

3.3.4  Transient  Waves 

Several  issues  arise  with  the  use  of  either  a  sinusoid  or  Gaussian  incident  wave.  The 
introduction  of  the  sinusoid  wave  at  the  surface  of  the  scatterer  produces  transient  waves  in  the 
computational  space.  The  frequency  data  for  the  RCS  calculations  cannot  be  taken  until  the 
transients  diminish.  The  electrical  size  of  the  target  determines  the  time  required  for  the 
transients  to  diminish.  Taflove  reports  that  the  minimum  time  is  four  times  the  longest  electrical 
length  of  the  object  [52].  The  required  time  may  be  shorter  or  longer,  as  will  be  seen  in  Chapter 
4,  to  obtain  converged  solutions. 
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3.3.5  Creeping  Waves  and  Traveling  Waves 

Creeping  waves  and  traveling  waves  are  electromagnetic  phenomena  that  can 
significantly  affect  the  RCS  results.  Knott  defines  creeping  waves  as  the  phenomenon  associated 
with  smooth  bodies  such  as  a  sphere  or  cylinder  which  are  not  large  with  respect  to  X.  The  wave 
is  launched  at  the  shadow  boundary  and  travels  around  on  the  shadow  side  of  the  object  and 
emerges  at  the  other  shadow  boundary  [17].  Because  the  amplitude  of  the  creeping  wave 
attenuates  exponentially,  creeping  waves  are  usually  negligible  for  bodies  which  are  larger  than 
IX  [27].  Since  the  cone-sphere  and  ogive  bodies  used  for  this  thesis  research  are  less  than  IX, 
creeping  waves  contribute  to  the  scattering  results  [17]. 

Surface  traveling  waves  are  an  electromagnetic  phenomenon  similar  to  creeping  waves 
except  the  wave  is  launched  along  smooth  surfaces  when  the  grazing  incidence  is  small.  The 
traveling  wave  occurs  on  the  “illuminated”  side  of  the  object,  as  opposed  to  the  creeping  wave 
which  occurs  on  the  shadow  side  of  an  object.  The  E  field  must  be  parallel  to  the  surface  and  a 
surface  discontinuity  must  exist  to  reflect  the  wave  back  along  the  surface. 

FVTD  will  accurately  calculate  creeping  and  surface  traveling  waves  if  the  grid  is 
correctly  generated  and  the  grid  density  at  the  surface  is  relatively  high.  The  waves  travel  close 
to  the  surface,  and  thus,  the  grid  must  be  dense  to  accurately  predict  the  propagation  of  waves 
near  the  surface.  The  exact  grid  density  required  will  be  discussed  in  more  detail  in  Chapter  4. 

3.3.6  Diffraction 

The  correct  prediction  of  diffraction,  another  electromagnetic  phenomenon,  is  critical  for 
the  accurate  calculation  of  scattering  results.  Diffraction  occurs  at  the  tips  or  edges  of  an  object. 
The  grid  requirements  for  correctly  predicting  diffraction  are  similar  to  the  grid  requirements  for 
creeping  and  traveling  waves.  The  requirements  for  the  grid  will  be  discussed  in  Chapter  4. 
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The  terms,  traveling  \vaves,  creeping  waves,  and  diffraction,  originated  with  the  use  of 
high  frequency  techniques  which  do  not  necessarily  account  for  these  interactions.  The  terms 
refer  to  electromagnetic  phenomena  that  have  to  be  accounted  for  by  using  other  methods  or 
correction  terms.  FVTD’s  direct  solution  of  Maxwell’s  equations  account  for  all  of  these 
phenomena  because  the  wave  physics  and  interactions  are  inherent  in  the  governing  equations. 
The  traditional  high  frequency  terms  are  used  in  the  discussion  to  easily  refer  to  the  phenomena. 

3.4  Comparisons/Benchmarks 

3.4.1  Method  of  Moments 

Analytical  solutions  do  not  exist  for  objects  such  as  the  ogive  and  cone-sphere.  A  low 
frequency  CEM  method  commonly  used  for  validation  purposes  is  the  Method  of  Moments 
(MoM).  The  accuracy  of  the  MoM  permits  it  to  be  used  for  validation.  The  MoM  requires  the 
surface  of  the  PEC  scatterer  to  be  discretized;  whereas,  for  FVTD,  a  volume  discretization  is 
required. 

The  MoM  solves  for  the  equivalent  current  density  induced  on  the  surface  of  the 
scatterer.  The  method  produces  a  system  of  linear  equations  which  can  be  written  in  matrix  form 
as 

V  =  Z-I  (22) 

where  I:  Current  Vector 

V :  V  oltage  V  ector 

Z:  Impedance  Matrix 

The  known  variables  are  the  voltage  vector  and  the  entries  of  the  impedance  matrix.  The  current 
vector  is  the  unknown;  therefore,  the  left-hand  inverse  of  Z  is  required  to  solve  for  I.  The  Z 
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matrix  is  usually  dense  and  can  become  very  large  for  electrically  large  objects,  limiting  the  size 
of  the  scatterer  that  can  be  studied.  The  inversion  of  the  Z  matrix  required  to  solve  for  the 
equivalent  currents  can  be  computationally  extensive  limiting  the  size  of  the  usable  matrix  [51]. 
The  traditional  radiation  integral  uses  the  equivalent  currents  to  calculate  the  scattered  fields. 

3.4.2  Empirical  Data 

The  EMCC  has  published  empirical  results  for  the  ogive  and  cone-sphere  discussed 
previously.  The  FVTD  computational  results  are  compared  to  the  monostatic  RCS  experimental 
data.  The  Volakis  reference  contains  the  details  on  the  measurement  methods  for  the  validation 
objects  [55]. 

3.4.3  Error  Calculations  and  Metrics 

The  FVTD  results  are  compared  to  MoM  and  empirical  RCS  data.  Metrics  are 
established  to  compare  the  FVTD  results  to  MoM  results  and  empirical  results.  Only  theoretical 
solutions  are  known  for  a  few  simple  shapes  such  as  the  sphere.  The  MoM  code  results  and  the 
empirical  data  for  the  ogive  and  cone-sphere  will  not  be  considered  “exact,”  but  as  very  accurate 
RCS  results  which  can  be  used  as  benchmarks. 

Four  methods  are  used  to  analytically  determine  the  difference  between  the  FVTD  and 
MoM  or  empirical  results.  The  various  methods  include  the 

•  difference,  in  dB,  after  considering  phase  by  calculating  the  minimum  error  between 
the  FVTD  RCS  result  at  one  angle  and  several  surrounding  (-I-/-5  degrees)  MoM  RCS 
results 

•  Mean-Square  Error  (MSE)  using  FVTD  and  MoM  results  (square  meters).  The  phase 
is  taken  into  account  in  the  same  way  as  for  the  previous  method 
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•  Cross-Correlation  between  FVTD  and  MoM  results  (square  meters) 

•  FFT  of  the  FVTD  and  MoM  results  (square  meters) 

The  first  two  methods  are  common  ways  to  calculate  the  difference  or  error.  As  will  be  seen  in 
several  plots  in  Chapter  4,  a  small  phase  shift  of  one  or  two  degrees  will  result  in  a  greater  error 
if  the  phase  shift  is  not  considered.  The  phase  shift  is  considered  for  the  first  two  methods  by 
calculating  the  minimum  difference  between  the  FVTD  RCS  result  at  one  point  and  the  MoM 
RCS  result  at  surrounding  points  (+/-5  degrees).  The  expression  for  the  MSE  is 

=  (23) 

n=l 

where  Xf:  FVTD  RCS  value  at  a  bistatic  or  monostatic  angle 

Xn,:  MoM  RCS  value  at  a  bistatic  or  monostatic  angle 

The  last  two  methods  are  mathematically  more  rigorous  techniques  to  compare  the 
FVTD  and  MoM  results.  The  cross-correlation  gives  the  similarity  between  two  plots  and 
locates  the  point  of  maximum  correlation.  The  absolute  value  of  the  Fourier  Transform 
eliminates  the  phase  shift  by  transforming  the  data  (in  degrees)  to  a  non-physical  quantity  in 
units  of  1/degrees.  A  MATLAB  program  which  calculates  the  difference  in  dB,  MSE,  and  the 
data  for  the  cross-correlation  and  Fourier  transform  plots  is  listed  in  Appendix  B. 

3.5  Computer  Support 

The  FVTD  code,  written  in  FORTRAN  77,  was  run  on  Air  Force  Institute  of  Technology 
(AFIT)  machines,  such  as  the  SUN  Sparc  20  and  Digital  Electric  Corporation  (DEC)  Alpha,  to 
obtain  familiarization  with  the  code.  After  modifications  were  made  to  the  code,  accounts  were 
created  on  high  performance  supercomputers  managed  by  high  performance  computing  centers 
(HPCC). 
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3.5.1  AFIT’s  Sparc  20 

AFIT’s  Sparc  20s  use  the  UNIX  operating  system.  The  Spares  were  useful  for  short 
FVTD  tests  because  of  the  relatively  quick  turn  around  time  and  accessibility.  For  tests  requiring 
longer  simulation  times  and  more  memory,  other  computing  machines  were  more  useful. 

3.5.2  AFIT’s  DEC  Alpha 

The  Digital  Electric  Corporation  (DEC)  Alpha  machines  use  an  AXP  central  processing 
unit  (CPU).  The  CPU  uses  a  21064  chip,  a  64-bit  reduced  instruction  set  chip  (RISC)  processor, 
which  operates  at  190  MHz.  The  DEC  Alpha  machines  utilize  the  OpenVMS  version  6.2 
operating  system.  The  DEC  Alpha  machines  were  useful  for  short  tests  for  the  3-D  objects. 

3.5.3  CEWES  HPCC’s  Cray  90 

The  Cray  90,  located  at  USAE  Waterways  Experiment  Station  (CEWES)  HPCC  in 
Vicksburg,  Mississippi,  is  a  16  processor  high  performance  vector  machine.  The  Cray  90,  or 
C90,  uses  the  UNICOS  9.0  operating  system.  The  peak  performance  for  each  processor  is  one 
Gigaflop.  The  FVTD  code,  optimized  for  a  vector  machine,  has  a  faster  ran  time  on  the  Cray  90 
than  any  other  computing  machine. 

3.5.4  Maui  HPCC’s  IBM  SP-2 

The  IBM  SP-2  is  a  scalable  parallel  machine  with  400  RS/6000  processors  which  operate 
at  66.7  MHz  each.  The  SP-2  is  a  high  performance  computing  machine  with  a  peak  throughput 
of  266  Mflops  per  processor.  The  IBM  SP-2  used  in  this  thesis  research  is  located  in  Maui, 
Hawaii,  at  the  Maui  HPC.  Although  the  FVTD  code  is  optimized  for  a  vector  machine,  shorter 
simulations  were  completed  on  one  processor  of  the  parallel  SP-2  machine. 
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3.6  Summary 

The  successful  code  modifications  permit  the  analysis  of  the  electromagnetic  scattering 
from  an  ogive  and  cone-sphere.  The  modified  code  can  be  applied  to  other  closed-surface 
perfect  electric  conductor  (PEC)  3-D  objects  if  an  appropriate  grid  is  generated.  In  addition  to 
modifying  the  code  for  a  generic  3-D  object,  the  option  of  specifying  the  direction  and 
polarization  of  the  incident  wave  was  programmed.  Other  options,  such  as  convergence  checks, 
ensure  accurate  RCS  data  is  computed. 

In  the  next  chapter,  the  options  programmed  into  the  code  are  used  to  perform  computer 
simulations  for  the  ogive  and  cone-sphere  for  different  grid  sizes,  frequencies,  types  of  incident 
waves,  etc.  The  option  added  to  specify  the  direction  and  polarization  of  the  incident  wave 
permits  the  completion  of  monostatic  tests.  The  monostatic  results  are  compared  to  MoM  and 
empirical  data.  The  different  computational  issues  discussed  in  this  chapter  serve  as  a 
foundation  for  explaining  the  results  and  the  differences  between  the  FVTD  and  MoM  or 
empirical  RCS  data.  The  ogive  and  cone-sphere  RCS  results  in  the  next  chapter  validate  the 
FVTD  code  and  algorithm  for  CEM. 
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4  Applications,  Results,  and  Comparisons 


4.1  Overview 

The  electromagnetic  scattering  and  RCS  results  for  the  EMCC-defined  RCS  test  bodies 
ogive  and  cone-sphere  are  presented  in  this  chapter.  Bistatic  and  raonostatic  RCS  results  are 
presented  for  each  PEC  test  body.  The  bistatic  results  are  compared  to  MoM  results  and  the 
monostatic  results  are  compared  to  empirical  data  and  MoM  data  to  validate  the  FVTD  algorithm 
and  code  for  CEM.  Sinusoid  incident  wave  results  for  the  ogive  are  compared  to  Gaussian  pulse 
incident  wave  results  to  illustrate  the  trade-offs  between  the  two  types  of  incident  waves.  The 
metrics  established  in  Chapter  3  depict  the  differences  between  the  FVTD  RCS  results  and  MoM 
and  empirical  RCS  data. 

To  analyze  the  scattering  from  the  ogive,  the  size  of  the  grid  is  varied  in  each  coordinate 
direction  (R,9,4))  to  perform  a  grid  convergent  study  and  to  obtain  the  optimum  grid  point  density 
(GPD)  for  each  coordinate  direction.  The  grid  for  the  cone-sphere  was  then  generated  using 
these  optimal  GPDs,  The  RCS  results  for  the  cone-sphere  confirm  the  grid  requirements 
obtained  for  the  ogive  and  validate  the  FVTD  algorithm  for  another  PEC  test  body. 

Several  different  computing  platforms  were  used  (See  Chapter  3)  to  complete  the  FVTD 
tests  for  each  scattering  object.  Several  test  simulations  were  completed  on  each  platform  to 
ensure  the  results  were  independent  of  the  platform.  Platform  differences  such  as  compiler 
optimization  and  floating  point  accuracy  can  affect  the  results.  For  example,  the  Cray  90  has  a 
floating  point  accuracy  of  14  digits  as  compared  to  the  7  digit  accuracy  of  the  DEC  Alpha  and 
SP-2  computing  machines.  Compiler  options  were  used  on  the  DEC  Alpha  and  SP-2  machines  to 
increase  the  floating  point  accuracy  to  that  comparable  to  the  Cray  90. 
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4.2  Ogive  Electromagnetic  Scattering  Results 


Table  1  is  the  test  matrix  listing  the  computer  simulations  completed  for  the  ogive.  The 
five  ogive  tests  include  several  subtests.  “OG”  in  the  test  number  refers  to  a  test  for  the  ogive, 
and  the  number  in  the  test  designator  refers  to  the  test  number.  The  last  letter  in  each  test 
number  refers  to  the  subtest.  The  subtests  are  groups  of  tests  that  use  a  specific  incident  wave, 
frequency,  or  grid  size.  In  the  discussion,  a  designator  such  as  OGIX  refers  to  the  entire  group 
of  subtests.  The  tests  are  easily  referenced  in  the  discussion  with  the  use  of  the  test  numbers. 


Table  1:  Ogive  Test  Matrix 


Test 

Incident  Wave 

Frequency 

Angle  of 

Grid  Size 

Number 

Type 

(GHz) 

Incidence 

OGla 

Sinusoid 

1.18 

0.0 

36-74-45 

OGlb 

(« 

(( 

71-43-45 

OGlc 

(« 

(( 

(( 

71-74-25 

OGld 

(« 

(i 

71-125-55 

OGle 

(( 

71-43-25 

OGlf' 

n 

(( 

71-74-45 

OG2a 

Gaussian 

1.18 

0.0 

71-74-45 

OG2b 

(« 

a 

OG3a‘ 

Sinusoid 

1.18 

10.0 

71-74-45 

OG3b' 

20.0 

(( 

OG3c‘ 

<( 

« 

30.0 

<( 

OG3d' 

40.0 

(C 

OG3e' 

n 

50.0 

(( 

(( 

60.0 

(C 

70.0 

(C 

OG3h' 

“ 

80.0 

(( 

OG3i‘ 

(( 

90.0 

(( 

OG4a 

Sinusoid 

9.0 

0.0 

61-125-95 

OG5a 

Gaussian 

1. 0-7.0,  Af=0.5 

0.0 

76-125-75 

OG5b 

Sinusoid 

1.0 

0.0 

83-74-35 

OG5c 

(( 

2.0 

« 

70-74-35 

OG5d 

(( 

3.0 

« 

69-74-45 

OG5e 

(( 

4.0 

tc 

76-125-75 

'Separate  tests  were  completed  for  HH  (transmit  horizontal,  receive  horizontal) 
and  VV  (transmit  vertical,  receive  vertical)  polarization.  The  grid  size  for  HH 
polarization  is  (71-74-45).  For  VV,  the  grid  size  is  (71-43-25). 
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The  RCS  results  for  the  ogive  for  each  test  is  compared  to  MoM  RCS  results  and 
experimental  data  for  VV  and  HH  polarization.  The  first  ogive  tests,  OGIX,  are  simulations 
which  analyze  the  variation  of  the  number  of  cells  in  each  coordinate  direction  (R,0,(t))  using  a 
sinusoid  incident  wave  at  1.18  GHz.  The  ogive  is  one  wavelength  long  at  1.18  GHz.  Tests 
OG2X  are  two  tests  at  1.18  GHz  using  a  Gaussian  pulse  for  the  incident  wave  with  the  frequency 
of  interest  near  the  center  of  the  pulse  and  near  the  edge  of  the  pulse.  Tests  OG3X  are  tests  for 
monostatic  calculations  at  1.18  GHz.  The  bistatic  RCS  for  the  ogive  at  9.0  GHz  is  test  OG4a. 
The  final  tests,  OG5X,  compare  results  using  sinusoid  incident  waves  for  various  frequencies  to 
the  results  from  a  Gaussian  pulse  incident  wave  test  to  illustrate  the  trade-offs  between  the  use  of 
each  incident  wave.  For  all  of  the  bistatic  tests,  the  angle  of  incidence  is  tip-on  at  0°. 

The  data  in  Table  2  lists  the  stability  and  time  data  for  each  computer  simulation  such  as 
the  CFL  value,  time  step,  number  of  time  steps  per  period,  total  number  of  periods,  total  number 
of  time  steps,  and  total  CPU  time.  The  CFL  value  used  for  the  majority  of  the  tests  is  1.5.  The 
maximum  CFL  value  to  enforce  stability  is  1.74;  however,  a  smaller  CFL  value  provides  greater 
accuracy.  A  value  of  1.5  provides  reasonable  accuracy  and  a  time  step  close  to  the  maximum. 
The  time  steps  are  normalized  to  the  wave  velocity  of  the  medium.  The  time  steps,  if  not 
normalized,  are  in  the  picosecond  range.  As  discussed  in  the  thesis  scope  in  Chapter  1,  this 
thesis  research  did  not  focus  on  speed  optimization.  Due  to  the  structure  of  the  grids  and  the  tips 
of  the  ogive,  the  cell  size  can  be  very  small  relative  to  X  and  the  resulting  time  step  is  very  small 
(i.e.  2.3221E-5  or  0.07740  picoseconds  for  test  OG4a).  The  total  number  of  periods  required  for 
the  test  is  related  to  the  electrical  size  of  the  object.  Normally,  the  computer  test  must  be,  in 
periods,  two  to  four  times  the  electrical  length  of  the  object.  The  computer  simulation  times  for 
several  tests  are  large  (e.g.  OG4a)  for  the  higher  frequencies  due  to  the  electrical  length  of  the 
object  and  the  small  time  step. 
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Table  2:  Stability  and  Time  Data  for  the  Ogive  Tests 


Test 

Number 

CFL 

At 

(Seconds(c)) 

Time  Steps/ 
Period 

Total 

Periods 

OGla 

1.5 

1.5030E-4 

1692 

7 

11844 

13194 

OGlb 

“ 

1.1674E-4 

2178 

« 

15246 

18137 

OGlc 

“ 

2.0210E-4 

1258 

8806 

8831 

OGld 

(« 

9.5292E-4 

2668 

it 

18676 

66744 

OGle 

(( 

2.0210E-4 

1258 

66 

8806 

N/A 

OGlf 

(( 

1.1674E-4 

2178 

66 

15246 

27773 

OG2a 

(( 

2.0210E-4 

N/A 

N/A 

14263 

N/A 

OG2b 

a 

2.0210E-4 

N/A 

N/A 

14263 

N/A 

OG3X-HH 

<t 

1.1674E-4 

2178 

7 

15246 

27851 

OG3X-VV 

1.7 

2.2904E-4 

979 

5 

4895 

3156 

1.7 

5.1749E-5 

644 

22 

14168 

79922 

OG5a 

1.5 

5.8723E-5 

N/A 

N/A 

17020 

117396 

OG5b 

(t 

1.4938E-4 

2009 

7 

14063 

22831 

OG5c" 

(( 

1.3159E-4 

1140 

8 

18240 

25902 

«( 

9.71272-5 

1030 

13 

13390 

24190 

OG5e" 

<( 

5.8723E-5 

1277 

14 

17878 

104909 

'Cray  90  CPU  time 

^he  RCS  convergence  check  was  used  for  these  FVTD  tests 


The  grid  point  densities  (GPD)  in  each  coordinate  direction  (R.O.tj))  for  each  ogive  test 
are  shown  in  Table  3.  The  grid  point  density  is  defined  as  the  number  of  finite-volume  cells  per 
wavelength,  X,  and  is  significant  when  generating  a  grid  for  an  object  to  ensure  accuracy.  In  the 
radial  direction,  as  explained  in  Chapter  3,  the  spacing  of  the  cells  increases  as  the  distance 
increases  from  the  surface  of  the  scatterer.  Therefore,  the  grid  point  density  is  much  larger  at  the 
surface  of  the  scatterer  as  compared  to  the  outer  edge  of  the  grid.  As  can  be  seen  in  the  table,  the 
GPD  varies  from  200  to  602  cells/A.  at  the  surface  and  from  8.2  to  24.7  cells/A,  at  the  outer  edge 
of  the  grid.  The  spacing  required  at  the  surface  is  due  to  the  object’s  small  electrical  size  and  the 
tips  which  requires  the  correct  calculation  of  the  electromagnetic  phenomena  occurring  at  the 
surface  such  as  creeping  waves,  traveling  waves,  and  diffraction.  In  the  radial  direction,  the 
outer  edge  of  the  grid  is  three  wavelengths  from  the  surface  of  the  scatterer  for  every  test. 
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Table  3:  Grid  Point  Densities  for  the  Ogive  Tests 


Test 

Number 

GPD;  Rat 
Surface 

GPD;  Rat 
Outer  Edge 

WilMM 

OGla 

255 

10.5 

57.4 

104 

64.2 

1975 

OGlb 

513 

21.1 

21.8 

104 

65.2 

1974 

OGlc 

510 

20.9 

57.4 

104 

32.1 

987 

OGld 

509 

20.8 

116 

104 

79.8 

2468 

OGle 

513 

21.1 

21.8 

104 

32.6 

987 

OGlf 

510 

20.9 

57.4 

104 

64.2 

1974 

OG2a 

513 

21.1 

21.8 

104 

32.6 

987 

OG2b 

513 

21.1 

21.8 

104 

32.6 

987 

OG3X-HH 

510 

20.9 

57.4 

104 

64.2 

1974 

OG3X-VV 

513 

21.1 

21.8 

104 

32.6 

987 

OG4a 

200 

8.2 

15.2 

13.6 

18.8 

582 

OG5a' 

375 

15.4 

34.2 

30.7 

32.9 

1019 

OG5b 

602 

24.7 

67.7 

123 

56.8 

1747 

OG5c 

451 

18.5 

33.8 

61.4 

28.4 

874 

OG5d 

401 

16.4 

22.6 

40.9 

25.2 

777 

OG5e 

375 

15.4 

34.2 

30.7 

32.9 

1019 

'Grid  optimized  for  4.0  GHz 


For  the  other  coordinate  directions,  theta  and  phi,  the  GPD  is  usually  the  smallest  at  the 
middle  of  the  ogive  and  largest  at  the  tips.  For  the  cases  in  which  the  theta  GPD  at  the  tips  is 
slightly  less  than  the  GPD  at  the  body,  as  for  test  OG5a,  the  original  grid  was  not  changed  in  the 
theta  direction  and  the  grid  spacing  at  the  tip  was  slightly  greater  in  the  middle  of  the  ogive  body 
as  compared  to  the  tip.  The  theta  GPD  at  the  tips  varies  from  13.6  to  123  cells/^.  In  the  phi 
direction,  the  body-of-revolution  object  forces  the  GPD  to  be  large  at  the  tips  and  smaller  in  the 
center  of  the  body.  The  phi  GPD  in  the  center  of  the  body  is  as  large  as  79.8  cells/?i  and  as  small 
as  18.8  cells/X.  The  GPD  at  the  tips  is  the  parameter  used  for  determining  the  optimal  GPD  in 


the  theta  direction  and  the  lowest  GPD  in  the  phi  direction  is  the  parameter  used  for  determining 
the  optimal  GPD  in  that  direction. 


4.2.1  Ogive  Bistatic  RCS  Results  for  1.18  GHz,  Sinusoid  Incident  Wave 

The  first  set  of  tests,  OGIX,  study  the  scattering  from  the  ogive  at  1.18  GHz  using  a 
sinusoid  incident  wave.  The  incident  wave  is  incident  at  0°,  or  tip-on.  Figures  12  through  15  are 
plots  of  the  nose-on  bistatic  scattering  from  the  ogive  for  VV  polarization.  Figures  12,  13,  and 
14  show  the  results  as  the  grid  point  density  is  varied  in  the  radial,  theta,  and  phi  directions, 
respectively  (i.e.  tests  OGla,  OGlb,  and  OGlc).  Test  OGlf  is  considered  the  acceptable  FVTD 
result  after  a  grid  convergence  study  was  completed.  The  original  size  of  the  grid  limited  the 
grid  sizes  which  could  be  generated.  The  number  of  grid  points  have  to  be  a  multiple  of  the 
original  grid  size. 

As  can  be  seen  in  Figure  12,  the  largest  errors  occur  when  the  number  of  cells  is 
decreased  in  the  radial  direction.  The  primary  reason  for  the  change  is  due  to  the  first-order 
surface  boundary  condition  which  does  not  accurately  predict  diffraction  and  traveling  waves 
unless  the  grid  points  are  tightly  packed  at  the  surface.  The  GPD  for  OGla  as  shown  in  Figure 
12  is  255  cells/X,  in  the  radial  direction  and  does  not  produce  good  results.  To  maintain  accuracy, 
the  first  two  cells  in  the  radial  direction  must  be  small  relative  to  the  wavelength.  For  the  ogive 
tests,  OGlb  through  OGlf,  the  grid  point  density  (GPD)  at  the  surface  in  the  radial  direction  is 
513  cells/^.  If  the  object  is  electrically  larger,  the  GPD  at  the  surface  in  the  radial  direction  does 
not  have  to  be  as  large,  as  will  be  seen  for  test  OG4a,  because  the  electromagnetic  phenomena 
can  be  considered  local. 

In  the  theta  and  phi  directions,  the  density  can  be  decreased  to  21.8-32.6  cells/?i  from  80- 
110  cells/A,  as  can  be  seen  in  Figure  15  to  obtain  accurate  results.  As  the  grid  point  density 
decreases,  the  errors  in  the  RCS  occur  first  in  the  backscatter  and  forward  scattering  directions  as 
would  be  expected  because  of  the  diffraction  at  the  tips  of  the  ogive. 
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In  each  plot  for  the  ogive  at  1.18  GHz,  VV,  the  FVTD  results  are  plotted  against  MoM 
results.  A  surface  grid  with  31  points  along  the  arc  was  used  to  obtain  the  MoM  results.  Figure 
16  shows  the  difference  or  error  between  the  FVTD  results  for  test  OGlf  and  the  CICERO  MoM 


Figure  14:  Ogive  Bistatic  RCS,  1.18  GHz,  W,  Number  of  Cells  Varied  in  (j)  Direction 


results.  The  difference  is  no  more  than  0.71  dB  if  the  phase  difference  is  considered  as  described 
in  Chapter  3.  The  largest  difference  occurs  in  the  forward  scattering  region.  Figure  17  illustrates 
the  cross-correlation  between  the  FVTD  results  from  test  OGlf,  VV  polarization,  and  MoM 
results.  The  plot  clearly  shows  the  FVTD  and  MoM  results  have  a  high  correlation  but  the  phase 
of  the  FVTD  result  is  one  or  two  degrees  to  the  left  of  the  MoM  results.  If  the  null  in  Figure  15 
occurred  at  the  same  bistatic  angle,  the  maximum  correlation  would  occur  at  180  degrees. 
Another  representation  of  the  accuracy  of  the  FVTD  results  is  shown  in  Figure  18  where 
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Difference  (dB)  ^  RCS  (dBsm) 

0.00  0.20  0.40  0.60  0.80  S  -60  -50  -40  -30 


e  15:  Ogive  Bistatic  RCS,  1.18  GHz,  W,  Fine  (71-125-55 )  vs.  Coarse  (71-43-25)  Grid 


Figure  16:  Ogive  Bistatic  RCS,  1.18  GHz,  W,  Difference  Between  FVTD  and  MoM  RCS 
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Theta  (Degrees) 

Figure  17:  Ogive  Bistatic  RCS,  1.18  GHz,  W,  Cross-Correlation  ofFVTD  and  MoM  RCS 

the  Fourier  transform  of  each  plot  is  taken.  The  x  axis  has  units  of  1/degrees  which  has  no 
physical  meaning.  The  maximum  magnitude  corresponds  to  the  large  140  degree  arc  in  the 
bistatic  plot.  Before  the  Fourier  transform  was  taken,  the  data  was  windowed  with  a  Hamming 
window.  If  the  FVTD  and  MoM  results  are  not  windowed  first,  the  truncation  of  the  data  in  the 
FT  sequence  produces  unwanted  oscillations  in  the  curve.  The  metrics  displayed  in  Figures  15- 
17  illustrate  the  small  difference  and  high  correlation  between  the  FVTD  and  MoM  RCS  results. 

The  frequency  data  for  tests  OGla  to  OGlf  was  taken  from  the  fifth  to  the  seventh 
periods  to  calculate  the  RCS  data.  Subsequent  tests  using  the  RCS  convergence  check  showed 
the  same  results  could  be  acquired  if  the  data  was  taken  from  the  fourth  to  the  fifth  periods. 
These  tests  revealed  that  for  1.18  GHz,  the  transients  introduced  with  the  sinusoid  incident  wave 
require  at  least  four  periods  to  diminish  before  frequency  data  can  be  taken  for  the  RCS 
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calculations.  Taflove  recommends  that  at  least  four  times  the  electrical  length,  in  periods,  are 
required  [51].  The  ogive  results  show  that  this  approximation  is  appropriate  for  this  frequency; 
however,  fewer  periods  can  be  used  for  higher  frequencies  to  obtain  accurate  data. 

Similar  results  were  obtained  for  HH  polarization.  The  same  grid  sizes  were  tested  and 
the  results  were  similar.  Figure  19  compares  the  coarse  grid  (OGle)  to  the  fine  grid  (OGld). 
Figure  20  is  the  plot  of  the  accepted  result  (OGlf)  and  the  MoM  results.  Again,  the  results  for 
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1 /Theta  (1 /Degrees) 


Figure  18:  Ogive  Bistatic  RCS,  1.18  GHz,  W,  Fourier  Transform  ofFVTD  and  MoM  RCS 


OGld  are  accurate  within  several  dB,  but  for  comparison,  the  slightly  more  accurate  results  from 
test  OGlf  are  used  to  compare  to  the  MoM  results. 

Figure  21  is  the  difference,  in  dB,  between  the  FVTD  and  MoM  results  for  HH 
polarization.  The  phase  shift  in  the  null  was  accounted  for  or  the  difference,  in  dB,  would  be 
much  greater  if  taken  point  by  point.  The  largest  difference  of  approximately  2.0  dB  occurs  at 
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151°.  The  null  in  the  FVTD  solution  at  that  location  is  not  as  deep  as  in  the  MoM  solution.  If 
the  nulls  are  not  considered,  the  greatest  difference  is  approximately  0.7  dB.  The  cross¬ 
correlation  was  computed  for  the  FVTD  and  MoM  RCS  results  and  the  FVTD  results  were 
highly  correlated  with  the  MoM  results,  except  the  shift  of  the  FVTD  data  was  to  the  right  by  one 
or  two  degrees.  This  shift  of  two  degrees  can  be  seen  visually  in  Figure  20.  The  Fourier 
transforms  of  the  curves  were  taken  and  were  very  similar  just  as  for  the  VV  RCS  results. 

The  FVTD  results  for  the  ogive  at  1.18  GHz  using  a  sinusoid  incident  wave  were 
excellent  when  compared  to  MoM  results.  The  results  illustrated  that  for  electrically  small 
objects,  the  FVTD  algorithm  correctly  computes  diffraction  and  traveling  waves.  The  small 
electrical  size  also  allowed  a  relatively  large  time  step  resulting  in  a  reasonable  simulation  time. 
For  test  OGle,  the  Cray  90  CPU  time  was  approximately  3200  seconds.  For  higher  frequencies, 
as  shown  in  Table  2,  the  computer  simulation  times  are  much  longer. 


Figure  19:  Ogive  Bistatic  RCS,  1.18  GHz,  HH,  Fine  (71-125-55)  V5.  Coarse  (71-43-25)  Grid 
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4.2.2  Ogive  Bistatic  RCS  Results  for  1.18  GHz,  Gaussian  Pulse  Incident  Wave 

The  accuracy  of  the  Gaussian  pulse  incident  wave  is  illustrated  in  Figures  22  and  23. 
One-third  of  the  bandwidth  of  the  pulse  is  considered  to  be  adequate  to  take  frequency  data.  This 
is  illustrated  in  the  plots  for  tests  OG2a  and  OG2b.  Test  OG2a  placed  the  frequency  of  interest, 
1.18  GHz,  near  the  center  of  the  pulse  bandwidth.  A  Gaussian  pulse  with  a  duration  of  0.017831, 
normalized  to  c,  and  a  radian  frequency  delta  of  12.35693  rad/s  was  used  to  place  the  frequency 
near  the  center  of  the  pulse.  The  test,  OG2b,  placed  the  frequency  near  the  edge  of  the  usable 
bandwidth  by  using  a  Gaussian  pulse  with  a  duration  of  0.100000  and  a  radian  frequency  delta  of 
2.471386  rad/s.  The  results  are  almost  identical.  Appendix  B  contains  detailed  data  for 
specifying  the  parameters  for  the  Gaussian  pulse  incident  wave. 

The  vertical  polarization  case  is  plotted  in  Figure  22.  The  Gaussian  pulse  RCS  results, 
OG2X,  are  compared  to  MoM  RCS  results.  The  horizontal  polarization  case  is  shown  in  Figure 
23.  The  Gaussian  pulse  RCS  results,  OG2X,  match  the  RCS  calculations  of  test  OGle  which 
also  used  the  coarse  grid  size  of  (71-43-25).  The  results  for  the  grid  size  of  (71-43-25)  are  nearly 
as  accurate  as  the  results  for  the  grid  size  of  (71-74-45)  to  permit  use  of  the  smaller  grid.  The 
savings  in  computer  time  also  justifies  the  use  of  the  smaller  grid. 

The  comparison  of  tests  OGIX  and  OG2X  illustrate  the  accurate  results  which  can  be 
obtained  with  either  a  Gaussian  pulse  or  sinusoid  incident  wave.  The  computer  mn  times  for  the 
Gaussian  pulse  tests,  OG2X,  are  greater  than  the  times  for  the  sinusoid  (approximately  9600 
seconds  vs.  3200  seconds).  This  factor  of  three  will  not  always  be  the  relationship  between  run 
times.  The  electrical  length  of  the  object  and  the  grid  size  will  affect  the  computer  run  times. 

The  sinusoid  incident  wave  is  clearly  the  correct  incident  wave  to  use  if  the  RCS  is  desired  for  a 
single  frequency.  The  Gaussian  pulse  is  used  when  the  simulation  time  is  less  than  the  total  time 
required  for  the  individual  sinusoid  tests.  Test  OG5X  looks  at  these  trade-offs  in  more  detail. 
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Figure  22:  Ogive  Bistatic  RCS,  1.18  GHz,  W,  Gaussian  Pulse  Incident  Wave 


Figure  23:  Ogive  Bistatic  RCS,  1.18  GHz,  HH,  Gaussian  Pulse  Incident  Wave 
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4.2.3  Ogive  Monostatic  RCS  Results  for  1.18  GHz 

In  addition  to  bistatic  results,  the  FVTD  code  can  also  obtain  monostatic  data.  Multiple 
simulations  have  to  be  completed  to  obtain  monostatic  data  for  one  frequency  as  compared  to  one 
test  for  bistatic  data.  One  simulation  produces  a  bistatic  plot  for  0°  to  180°.  The  simulation 
produces  a  monostatic  result  for  one  angle,  the  angle  of  incidence.  To  obtain  a  full  monostatic 
sweep,  a  bistatic-to-monostatic  approximation  is  used  (See  Chapter  3).  The  approximation 
requires  tests  to  be  completed  every  10°  and  bistatic  data  completes  the  monostatic 
approximation.  Tests  OG3X  are  monostatic  calculations  for  the  ogive  at  1.18  GHz.  A  test  was 
completed  for  an  angle  of  incidence  every  10°  from  0°  to  90°.  The  ogive  is  symmetric  about  the 
xy  plane  resulting  in  a  symmetric  monostatic  plot  about  6=90°. 

The  monostatic  approximation  for  HH  polarization  is  plotted  in  Figure  24.  The  HH 
monostatic  test  used  a  grid  size  of  (71-74-45),  and  frequency  data  was  taken  from  the  fifth  to  the 
seventh  period.  The  FVTD  results  are  plotted  against  MoM  RCS  results  and  empirical  RCS  data. 
As  can  be  seen  in  the  plot,  the  MoM  and  FVTD  results  are  almost  identical  and  differ  from  the 
empirical  data  by  nearly  the  same  value.  In  Figure  25,  the  difference  between  the  FVTD,  MoM, 
and  empirical  results  is  plotted.  The  FVTD  results  differ  from  the  MoM  results  by  no  more  than 
2.5  dB.  If  the  large  fluctuations  are  ignored  in  the  empirical  data,  the  FVTD  results  are  within 
3.1  dB  of  the  empirical  data.  The  FVTD  results  differ  from  the  empirical  data  by  2.8  dB  at  the 
tips  (0=0°  and  0=180°)  and  2.5  dB  at  broadside  incidence  (6=90°). 

The  VV  polarization  results  for  the  monostatic  calculations  are  shown  in  Figure  26.  For 
this  group  of  tests,  a  grid  size  of  (71-43-25)  was  used.  Again,  the  FVTD  results  are  excellent 
when  plotted  against  MoM  and  experimental  RCS  data.  The  FVTD  results,  however,  do  not 
agree  well  at  the  nulls  at  30°  and  150°.  The  difference  between  the  FVTD  results  and  the  other 
data,  shown  in  the  figure,  can  be  as  large  as  8-10  dB  in  the  nulls. 
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The  accuracy  and  limitations  of  the  bistatic-to-monostatic  approximation  can  be 
determined  by  close  observation  of  the  FVTD  monostatic  plots.  The  data  at  10°  increments,  such 
as  0°,  10°,  20°,  etc.,  are  the  tme  monostatic  FVTD  RCS  results  and  are  in  agreement  with  the 
MoM  RCS  results.  However,  the  data  at  the  angles  5°,15°,25°,35°,45°,  etc.  are  at  the  junction  of 
the  approximation  data  from  each  monostatic  test.  Slight  discontinuities  can  be  seen  at  these 
angles.  For  example,  the  monostatic  data  for  25°-35°  degrees  is  obtained  from  the  bistatic  test 


Figure  26:  Ogive  Monostatic  RCS,  1.18  GHz,  W 


at  30°.  The  monostatic  data  for  35°-45°  is  taken  from  the  40°  bistatic  test.  A  small  discontinuity 
can  be  seen  at  35°.  The  approximation  is  only  good  for  small  angles,  5°  or  less.  If  the 
discontinuity  is  great,  the  approximation  is  not  poor,  but  more  monostatic  tests  have  to  be  made 
for  more  angles  so  the  approximation  is  used  for  smaller  angles.  A  large  number  of  oscillations 
in  the  plot  can  require  a  greater  number  of  monostatic  tests  (i.e.  higher  frequencies). 
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4.2.4  Ogive  Bistatic  RCS  Results  for  9.0  GHz 

The  ogive  tests  completed  previously  were  at  one  frequency,  1.18  GHz.  Further  analysis 
is  completed  for  the  ogive  using  FVTD  by  running  a  test,  OG4a,  for  the  ogive  at  9.0  GHz.  The 
ogive  is  approximately  7.6  wavelengths  long  at  this  frequency.  The  grid  size  required  for  this 
frequency  is  much  larger  along  with  an  increased  simulation  time.  The  wavelength  is  smaller  at 
9.0  GHz,  but  the  time  step  is  much  smaller  due  to  the  fine  grid  resulting  in  a  large  simulation 
time  (i.e.  79,922  seconds)  as  seen  in  Table  2. 

The  ogive  results  for  9.0  GHz,  HH  polarization,  are  shown  in  Figure  27.  The  FVTD 
results  for  9.0  GHz,  VV,  are  shown  in  Figure  28.  The  FVTD  results  are  plotted  against  MoM 
RCS  data.  The  results  are  excellent,  except  there  are  small  discrepancies  for  the  backscatter  and 
forward  scattering  regions.  The  grid  densities  for  these  tests  are  smaller  than  for  1.18  GHz.  The 
GPDs  in  the  theta  and  phi  directions  are  15.2  cells/X  and  18.8  cells/^,  respectively.  The  results 
for  the  ogive  at  1.18  GHz  showed  that  a  GPD  of  approximately  22-32  cells/A.  gives  the  best  data. 
These  results  depict  the  dependence  of  the  required  GPD  on  the  electrical  size  of  the  object.  The 
GPD  can  be  smaller,  15-20  cells/^,  as  the  electrical  length  of  the  object  gets  larger. 

The  9.0  GHz  results  for  the  ogive  illustrate  the  dependence  of  the  length  of  simulation 
time,  in  periods,  to  the  length  of  the  object.  At  1.18  GHz,  the  test  had  to  be  at  least  four  times,  in 
periods,  the  length  of  the  object.  The  same  factor  would  require  a  simulation  time  of  30  periods 
for  9.0  GHz.  This  is  not  required  because,  at  9.0  GHz,  the  ogive  is  in  the  optical  region.  As  a 
rule  of  thumb,  the  optical  region  is  encountered  for  an  object  greater  than  approximately  1%  [27]. 
At  9.0  GHz,  the  diffraction  and  the  traveling  waves  can  be  considered  to  be  more  of  a  local 
phenomena  than  for  1.18  GHz.  This  reduces  the  simulation  time  for  the  test  to  approximately 
three  times  the  length  of  the  object,  in  periods,  instead  of  four. 
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The  results  for  the  ogive  at  9.0  GHz  begin  to  show  the  limitations  of  the  explicit  FVTD 
algorithm  and  available  computing  resources.  The  stability  of  the  explicit  time-domain  algorithm 
is  related  to  the  cell  size  which  is  small  for  this  grid  at  9.0  GHz.  The  grid  size  is  large  to  meet 
the  GPD  requirements  and  the  fineness  of  the  grid  at  the  tips  forces  the  time  step  to  be  small. 

The  small  time  step,  the  large  grid  size,  and  the  number  of  periods  required  to  obtain  accurate 
frequency  data  results  in  a  long  simulation  time  (i.e.  79,922  seconds).  Accuracy  in  the 
backscatter  and  forward  scattering  regions  is  sacrificed  for  a  faster  run  time  and  smaller  grid. 

4.2.5  Ogive  Bistatic  RCS  Results  for  the  Gaussian  Pulse  Incident  Wave 

The  option  of  specifying  either  a  sinusoid  or  Gaussian  pulse  for  the  incident  wave  gives 
the  user  flexibility  in  obtaining  FVTD  RCS  data  for  one  frequency  or  multiple  frequencies, 
respectively,  depending  on  the  results  desired.  The  Gaussian  pulse  provides  RCS  data  for 
multiple  frequencies  but  normally  requires  a  longer  run  time  (See  Table  2);  therefore,  a  sinusoid 
is  ideal  for  single  frequency  RCS  data.  Included  in  this  section  are  results  for  several  frequencies 
using  a  Gaussian  pulse.  The  Gaussian  pulse  RCS  data  are  compared  to  sinusoid  incident  wave 
RCS  data  and  the  results  for  one  frequency,  4.0  GHz,  are  also  compared  to  MoM  RCS  results. 

The  Gaussian  pulse  used  for  test  OG5a  was  optimized  for  4.0  GHz.  The  grid  used,  (76- 
125-75),  was  also  optimized  for  that  frequency.  The  bistatic  RCS  results  for  frequencies  less 
than  4.0  GHz  are  excellent  and  will  be  shown.  The  GPDs  due  to  the  grid  size  for  frequencies 
below  4.0  GHz  are  sufficient  to  provide  excellent  results.  The  Gaussian  pulse  for  test  OG5a  had 
a  width  of  0.0182556  (See  Appendix  B)  with  a  radian  frequency  delta  of  10.471975  rad/s  (0.5 
GHz).  The  Gaussian  pulse  width  had  a  bandwidth  of  21.0  GHz  and  a  usable  BW  of  7.0  GHz  for 
these  parameters.  The  grid  was  optimized  for  4.0  GHZ  so  results  for  higher  frequencies  are  not 
discussed. 
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The  FVTD  results  for  1.0  GHz  using  the  Gaussian  pulse  (OG5a)  and  a  sinusoid  (OG5b) 
are  shown  in  Figure  29,  HH,  and  Figure  30,  W.  The  RCS  convergence  check  was  not  used  for 
test  OG5b.  The  results  are  identical  except  for  a  negligible  small  difference  in  the  forward 
scattering  region  which  could  be  due  to  a  different  grid  size,  numerical  dispersion,  or  reflections 
from  the  outer  boundary.  Slight  differences  were  seen  previously  for  different  grid  sizes.  The 
longer  run  time  for  the  Gaussian  pulse  test  can  introduce  numerical  dispersion  which 
accumulates  during  the  test.  Also,  the  longer  simulation  for  the  Gaussian  pulse  may  allow  the 
reflections  from  the  outer  boundary  to  affect  the  results  more  than  for  the  sinusoid  test. 

The  Gaussian  pulse  incident  wave  also  provided  RCS  data  for  2.0  GHz.  The  HH 
polarization  results  are  shown  in  Figure  31.  As  can  be  seen  in  the  figure,  the  data  for  the 
Gaussian  pulse  (OG5a)  and  sinusoid  incident  wave  (OG5c)  are  almost  identical.  A  small 
difference  exists  between  the  results  in  the  forward  scatter  direction.  The  RCS  convergence 
check  was  used  for  test  OG5c.  The  VV  polarization  result  is  shown  in  Figure  32.  The  results  are 
identical  except  for  a  slight  difference  in  the  forward  scatter  region. 

The  results  comparing  Gaussian  pulse  RCS  data  to  sinusoid  data  reveal  several  important 
computational  issues  requiring  attention  when  obtaining  FVTD  results.  The  Gaussian  pulse 
results  accurately  predict  the  backscatter  and  forwMd  scatter  regions.  Results  from  a  sinusoid 
incident  wave  will  not  if  the  transients  have  not  diminished.  The  convergence  check  for  the 
sinusoid  was  programmed  to  take  data  for  one  period,  sample  the  RCS  values,  take  data  for 
another  period,  and  sample  again.  Other  FVTD  tests  showed  that  more  accurate  results  can  be 
obtained  from  sampling  data  over  a  two  to  three  period  span.  A  convergence  check  designed  to 
sample  every  two  or  three  periods  would  provide  a  more  robust  and  accurate  method  of  ensuring 
accurate  results  from  a  sinusoid  incident  wave,  especially  for  higher  frequencies.  At  higher 
frequencies,  the  enormous  run  times  require  methods  which  ensure  accurate  results  with  one  test. 
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'igure  29:  Ogive  Bistatic  RCS,  1.0  GHz,  HH,  Gaussian  Pulse  vs.  Sinusoid  Incident  Wave 


Figure  30:  Ogive  Bistatic  RCS,  1.0  GHz,  W,  Gaussian  Pulse  V5.  Sinusoid  Incident  Wave 
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Figure  31:  Ogive  Bistatic  RCS,  2.0  GHz,  HH,  Gaussian  Pulse  vs.  Sinusoid  Incident  Wave 
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Figure  32:  Ogive  Bistatic  RCS,  2.0  GHz,  W,  Gaussian  Pulse  vs.  Sinusoid  Incident  Wave 
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Another  computational  issue  which  arises  with  the  use  of  the  Gaussian  pulse  is  the 
location  of  the  outer  boundary  for  each  frequency.  The  grid,  when  using  a  Gaussian  pulse,  is 
optimized  for  a  specific  frequency  with  the  location  of  the  outer  boundary  three  wavelengths 
from  the  scatterer.  The  test  gives  results  for  lower  frequencies,  but  the  outer  boundary  for  the 
lower  frequencies  is  less  than  three  wavelengths  away  from  the  scatterer  surface.  For  example, 
the  outer  edge  of  the  grid  for  2.0  GHz  in  test  OG5a  is  1.5  wavelengths  from  the  surface  of  the 
scatterer.  As  the  frequency  decreases,  the  reflections  from  the  outer  boundary  will  increase  the 
errors  for  those  particular  frequencies.  The  slight  errors  in  the  forward  scatter  regions  for  1.0 
and  2.0  GHz  are  most  likely  due  to  the  reflections  from  the  outer  boundary  as  opposed  to 
numerical  dispersion  in  the  Gaussian  pulse  test. 

In  Figure  33,  the  accurate  results  for  3.0  GHz,  HH,  is  illustrated  using  a  Gaussian  pulse 
(OG5a)  and  a  sinusoid  (OG5d).  The  VV  polarization  case  is  shown  in  Figure  34.  The  grid  sizes 
for  each  test  are  different,  but  the  results  are  almost  identical.  No  errors  in  the  forward  scatter 
region  reveal  that  numerical  reflections  are  causing  the  slight  errors  in  the  forward  scatter  area 
for  the  lower  frequencies  of  1 .0  and  2.0  GHz. 

Test  OGld  ran  for  thirteen  periods,  although  convergence  was  almost  reached  at  10 
periods  (within  0.12  dB  as  compared  to  0.1  dB  for  the  RCS  convergence  check).  At  10  periods, 
the  simulation  time  is  approximately  four  times  the  electrical  length  of  the  ogive.  The 
electromagnetic  phenomena,  including  diffraction  at  the  tips  and  traveling  waves  along  the 
surface  of  the  body,  require  time  to  stabilize.  These  frequencies  dictate  that  the  electromagnetic 
scattering  cannot  be  considered  local  phenomena.  Multiple  diffractions  from  the  tips  due  to 
traveling  waves  require  time  to  establish  and  stabilize.  As  the  electrical  size  of  the  object  gets 
larger,  as  for  the  case  of  9.0  GHz,  the  multiple  diffractions  are  not  as  prominent  and  less  time  is 
required  to  establish  the  stabilized  fields. 
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'igure  33:  Ogive  Bistatic  RCS,  3.0  GHz,  HH,  Gaussian  Pulse  vs.  Sinusoid  Incident  Wave 
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Figure  34:  Ogive  Bistatic  RCS,  3.0  GHz,  W,  Gaussian  Pulse  vs.  Sinusoid  Incident  Wave 


To  conclude,  the  excellent  results  using  a  Gaussian  pulse  at  4.0  GHz  for  the  ogive  is 
discussed.  Tests  OG5a  and  OG5e  were  optimized  for  4.0  GHz.  The  ogive  is  approximately  3.4 
wavelengths  long  at  this  frequency.  The  FVTD  RCS  data  using  the  Gaussian  pulse  are  compared 
to  MoM  RCS  results  in  Figure  35.  The  negligible  differences  are  in  the  nulls  and  in  the  forward 
and  backscatter  directions.  The  greatest  difference  between  the  FVTD  and  MoM  RCS  results  is 
approximately  1.5  dB  as  shown  in  Figure  36.  To  obtain  the  excellent  results,  a  GPD  for  theta 
and  phi  are  specified  at  30.7  cells/A,  and  32.9  cells/X,,  respectively.  Figure  37  shows  the  high 
correlation  between  the  MoM  and  FVTD  results.  The  nulls  in  the  bistatic  RCS  results  occur  at 
the  exact  same  angle.  The  cross-correlation  shows  the  highest  correlation  at  180°,  signifying  the 
best  correlation  possible. 

The  VV  polarization  results  for  the  ogive  at  4.0  GHz  are  shown  in  Figure  38.  The  FVTD 
results  are  plotted  against  MoM  results  and  show  excellent  agreement.  The  largest  difference  of 
1 .7  dB  occurs  at  the  monostatic  angle  of  0°  for  test  OG5a.  For  the  VV  results,  there  is  no  phase 
shift  which  would  result  in  a  high  cross-correlation  just  as  for  the  HH  results. 

To  complete  the  bistatic  and  monostatic  results  for  the  ogive,  an  example  of  the 
visualization  of  the  results  achievable  from  the  FVTD  time-domain  code  is  shown  in  Figure  39. 
The  plot  is  a  contour  slice  of  the  scattered  electric  field  from  the  ogive,  for  4.0  GHz,  in  the  xy 
plane.  The  plot  is  a  gray-shade  contour  plot  of  the  bistatic  field  with  the  angle  of  incidence  at  0° 
(+z  direction).  The  monostatic  angle  is  0°  and  the  forward  scatter  angle  is  180°  (-z).  The  darkest 
areas  represent  the  strongest  field  and  the  lightest  areas  correspond  to  the  weakest  fields.  The 
strongest  field  in  the  contour  plot  corresponds  to  the  peak  in  Figure  38  at  approximately  140 
degrees.  As  mentioned  previously,  the  ogive  is  3.4  wavelengths  long  at  this  frequency.  The 
same  number  of  wavelengths  can  be  seen  in  the  contours  along  the  edge  of  the  ogive. 
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'igure  35:  Ogive  Bistatic  RCS,  4.0  GHz,  HH,  Gaussian  Pulse  V5.  Sinusoid  Incident  Wave 
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Figure  36:  Ogive  Bistatic  RCS,  4.0  GHz,  HH,  Difference  Between  RCS  Results 
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Figure  37:  Ogive  Bistatic  RCS,  4.0  GHz,  HH,  Cross-Correlation  ofFVTD  and  MoM  RCS 


Several  important  observations  can  be  obtained  from  the  contour  plot  of  the  electric  field 
in  Figure  39.  The  perfect  symmetry  of  the  contours  about  the  z  axis  validates  the  computational 
accuracy  for  the  scattered  electric  field.  The  fields  in  each  half  of  the  plot  about  the  z  axis  are 
calculated  separately.  Any  errors  in  the  code,  especially  overlap  regions  in  the  grid,  would  create 
errors  in  the  symmetry  of  the  contour  plot.  The  symmetry  confirms  the  accuracy  of  the  code. 

The  diffraction  of  the  electromagnetic  wave  can  also  be  seen  at  the  tips  of  the  ogive.  In 
the  backscatter  region,  0°,  a  strong  field  exists.  The  grid  must  be  accurately  generated  for  the 
FVTD  algorithm  to  correctly  predict  the  diffraction  at  the  tip.  The  incident  field  and  the 
scattered  electric  and  magnetic  fields  are  implemented  at  cell  centers  and  not  cell  faces  or  cell 
grid  points.  The  grid  points  of  the  cells  define  the  tip  of  the  ogive,  but  if  only  the  cell  centers  are 
viewed,  the  tip  of  the  ogive  appears  to  have  a  small  blunt  tip.  This  characteristic  of  the  grid  and 
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Figure  38:  Ogive  Bistatic  RCS,  4.0  GHz,  VV,  Gaussian  Pulse  vs.  Sinusoid  Incident  Wave 
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Figure  39:  Ogive  Bistatic  Results,  4.0  GHz,  Electric  Field  Contour  Plot 


the  FVTD  formulation  demands  that  the  cells  at  the  tips  of  the  ogive  be  extremely  small  to 
accurately  represent  the  geometry  and  to  accurately  predict  the  diffracted  fields. 


To  finish  the  ogive  portion  of  the  validation,  two  parameters  are  calculated  to 
analytically  represent  the  excellent  RCS  results  for  the  scatterer  using  FVTD.  The  largest 
difference,  in  dB,  and  the  mean-square  error  (MSE),  explained  in  Chapter  3,  provide  metrics  for 
establishing  the  difference  and  error  between  FVTD  and  MoM  RCS  data.  Table  4  lists  the 
largest  dB  difference  and  MSE  for  each  polarization  of  each  ogive  test. 


Table  4:  MSE  and  Comparisons  for  the  Ogive  Tests 


Test  Number 

Difference, 
HH  (dB) 

MSE,  HH 

Difference, 
VV  (dB) 

MSE,VV 

OGld 

2.0 

1.80E-9 

0.8 

5.90E-11 

OGle 

3.0 

1.79E-9 

0.9 

5.92E-10 

OGlf 

2.0 

1.64E-9 

0.7 

7.38E-11 

OG2a 

3.0 

1.82E-9 

1.1 

4.50E-10 

OG2b 

2.9 

1.76E-9 

1.1 

4.34E-10 

OG3X 

2.5 

N/A 

N/A 

N/A 

OG4a 

2.4 

1.08E-5 

1.9 

8.52E-6 

OG5a  (4  GHz) 

1.4 

7.21E-8 

1.1 

4.74E-8 

OG5e 

1.4 

7.86E-8 

1.0 

5.30E-8 

The  largest  difference  between  bistatic  and  monostatic  FVTD  and  MoM  RCS  results  for 
the  tests  is  3.0  dB,  HH,  and  1.9  dB,  VV.  In  general,  the  errors  for  VV  are  less  than  the  errors  for 
HH  due  to  the  additional  nulls  seen  in  the  HH  plots.  The  results  showed  that  approximately  22- 
25  cells/A,  were  required  to  obtain  accurate  results  for  the  electrically  small  ogive.  Test  OG4a, 
the  9.0  GHz  test,  had  a  GPD  of  only  15.2  cells/A.,  illustrating  the  smaller  GPD  requirement  for 
electrically  larger  objects.  Test  OG4a  had  the  largest  MSE  and  test  OGlf,  HH,  had  the  smallest 
MSE.  These  MSE  results  are  visually  confirmed  in  the  RCS  plots.  The  small  MSEs  for  all  of 
the  tests  depict  the  excellent  results  obtained  for  the  ogive  using  FVTD. 
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4.3  Cone-Sphere  Electromagnetic  Scattering  Results 

The  RCS  calculations  for  the  cone-sphere  provide  further  validation  of  the  FVTD  code 
and  algorithm.  The  cone-sphere  is  a  common  RCS  test  body  but  the  narrow  cone  portion  and  the 
sphere  cap  provide  a  unique  body  for  analysis.  At  smaller  frequencies,  the  scattering  from  the 
cone-sphere  is  a  result  of  traveling  waves  along  the  narrow  cone,  creeping  waves  around  the 
sphere  cap,  and  diffraction  from  the  tip.  Excellent  RCS  results  for  the  cone-sphere  validate  the 
FVTD  code  and  algorithm  for  predicting  complicated  electromagnetic  phenomena. 

The  test  matrix  for  the  cone-sphere  computer  simulations  is  shown  in  Table  5.  Five  tests 
were  completed  which  include  several  subtests.  The  test  numbers  are  constructed  in  the  same 
method  as  the  ogive  test  numbers.  “CS”  is  used  instead  of  “OG”  to  refer  to  the  cone-sphere.  The 
test  numbers  are  used  in  the  discussion  of  the  results  to  easily  reference  the  tests. 

The  electromagnetic  scattering  results  for  the  cone-sphere  are  compared  to  MoM  RCS 
results  and  experimental  data  for  VV  and  HH  polarization.  The  first  cone-sphere  tests,  CSIX, 
are  bistatic  RCS  results  which  analyze  several  grid  sizes  at  0.869  GHz  using  a  sinusoid  incident 
wave  incident  on  the  sphere  section  of  the  cone-sphere.  The  cone-sphere  is  two  wavelengths 
long  at  0.869  GHz.  An  incident  angle  of  0°  corresponds  to  incidence  on  the  cone-sphere  along 
the  axis  of  symmetry  directly  onto  the  sphere-cap.  Tests  CS2X  are  bistatic  tests  at  0.869  GHz 
with  a  sinusoid  incident  wave  at  180°  or  tip-on  incidence.  Tests  CS3X  are  tests  for  the 
monostatic  calculations  at  0.869  GHz.  The  monostatic  RCS  is  computed  every  10°  from  0°  to 
180°  and  bistatic  data  completes  the  monostatic  plot.  The  cone-sphere  is  not  symmetric  about 
the  xy  plane  and  therefore  monostatic  tests  cannot  be  completed  for  only  0°  to  90°  as  for  the 
ogive.  CS4a  and  CS5a  are  bistatic  simulations  for  the  cone-sphere  at  3.0  GHz  using  a  sinusoid 
incident  wave  at  sphere-cap  (0=0°)  and  tip-on  incidence  (0=180°),  respectively. 
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Table  5:  Cone-Sphere  Test  Matrix 


Test  Number 

Frequency 

(GHz) 

Angle  of 
Incidence 

Grid  Size 

CSla' 

Sinusoid 

0.869 

0.0 

50-73-45 

CSlb' 

(( 

(( 

50-39-35 

CSlc 

« 

(( 

50-73-45 

CSld 

“ 

50-39-35 

CS2a* 

Sinusoid 

0.869 

180.0 

50-39-35 

CS2b 

“ 

a 

50-39-35 

CS3a 

Sinusoid 

0.869 

10.0 

50-39-35 

CS3b 

<( 

20.0 

a 

CS3c 

<( 

(( 

30.0 

n 

CS3d 

«( 

« 

40.0 

66 

CS3e 

<« 

50.0 

66 

CS3f 

u 

a 

60.0 

66 

CS3g 

<« 

(( 

70.0 

66 

CS3h 

(( 

80.0 

66 

CS3i 

n 

90.0 

(( 

CS3j 

« 

(( 

100.0 

(( 

CS3k 

110.0 

(( 

CS31 

«« 

(( 

120.0 

(( 

CS3m 

130.0 

(( 

CS3n 

140.0 

(( 

CS3o 

150.0 

(( 

CS3p 

160.0 

(( 

CS3q 

170.0 

(« 

CS4a‘ 

3.0 

0.0 

73-141-81 

CS5a' 

3.0 

180.0 

73-141-65 

'The  grid  spacing  at  the  surface  (R  direction)  for  these  tests  is  approximately  200  cells/A<. 
The  surface  grid  spacing  for  all  other  tests  is  400  cells/A,. 


The  data  in  Table  6  contains  the  stability  and  time  data  for  each  cone-sphere  test  such  as 
the  CFL  value,  time  step,  number  of  time  steps  per  period,  total  number  of  periods,  total  number 
of  time  steps,  and  total  CPU  time.  The  grid  is  structured  to  produce  the  largest  time  step  as 
possible.  For  example,  the  tip  GPD  in  the  theta  direction  is  22.0  cells/?i.  Greater  accuracy  could 
be  obtained  with  a  larger  GPD,  but  the  trade-off  between  run  time  and  accuracy  justifies  the 
spacing.  As  will  be  seen  in  the  discussion,  the  results  are  excellent  despite  the  larger  spacing  at 
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the  tip  of  the  cone-sphere.  Due  to  the  larger  time  step  as  compared  to  the  ogive,  the  majority  of 
the  tests  for  the  cone-sphere  were  completed  on  the  DEC  Alpha  and  SP-2  computing  machines. 


Table  6:  Stability  and  Time  Data  for  the  Cone-Sphere  Tests 


Test 

Number 

CFL 

At 

(Seconds(c)) 

Time  Steps/ 
Period 

Total 

Periods 

Total  Time 
Steps 

CPU  Time^ 
(seconds) 

m^m 

1.5 

1.9553E-4 

1766 

4 

7064 

9797 

n 

3.9937E-4 

864 

4 

3456 

2246 

(t 

1.7304E-4 

3990 

4 

15960 

22900 

ct 

2.9814E-4 

1158 

4 

4632 

3010 

CS2a^ 

n 

3.9937E-4 

864 

4 

3456 

2248 

CS2b^ 

a 

2.9814E-4 

1158 

4 

4632 

3072 

CS3X^ 

(( 

3.3789E-2 

1158 

4^ 

4632 

2697 

(( 

4.8238E-5 

2073 

9 

18657 

126612 

CS5a" 

(( 

5.3004E-5 

1887 

10 

18870 

102190 

'Cray  90  CPU  time 

^The  RCS  convergence  check  was  used  for  these  FVTD  simulations 
^Several  monostatic  tests  required  five  periods  to  reach  convergence 


Table  7:  Grid  Point  Densities  for  the  Cone-Sphere  Tests 


Test 

Number 

GPD:  Rat 
Surface 

GPD:  Rat 
Outer  Edge 

li 

GPD:  0  at 
Tip 

CSla 

201 

3.67 

51.8 

45.3 

29.6 

2367 

CSlb 

200 

3.65 

26.0 

22.0 

22.0 

862 

CSlc 

402 

2.64 

51.8 

45.3 

29.6 

2367 

CSld 

401 

2.56 

25.9 

22.0 

22.0 

862 

CS2a 

200 

3.65 

26.0 

22.0 

22.0 

862 

CS2b 

401 

2.56 

25.9 

22.0 

22.0 

862 

CS3X 

401 

2.56 

25.9 

22.0 

22.0 

862 

CS4a 

233 

9.3 

30.0 

26.2 

16.3 

2605 

CS5a 

233 

9.3 

30.0 

26.2 

14.0 

2056 

Table  7  shows  the  grid  point  densities  (GPD)  in  each  direction  for  each  cone-sphere  test. 
Again,  the  GPD  is  defined  as  the  number  of  finite-volume  cells  per  wavelength,  X,  and  must  be 
considered  when  generating  a  grid  for  a  scatterer.  As  explained  in  Chapter  3,  the  spacing  of  the 
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cells  in  the  radial  direction  increases  as  the  distance  from  the  surface  of  the  scatterer  increases. 
Therefore,  the  grid  point  density  is  much  larger  at  the  surface  as  compared  to  the  outer  edge  of 
the  grid.  As  can  be  seen  in  the  table,  the  GPD  varies  from  200  to  401  cells/A,  at  the  surface  and 
from  2.56  to  9.3  cells/A,  at  the  outer  edge  of  the  grid  for  the  cone-sphere  grid.  The  outer  edge  of 
the  grid  is  three  wavelengths  from  the  surface  of  the  scatterer  for  every  simulation.  The  GPD  in 
the  0  direction  is  approximately  the  same  at  the  sphere/cone  junction  and  the  cone  tip.  The  GPD 
for  the  tests  in  the  (j)  direction  is  the  smallest  at  the  sphere/cone  junction  and  the  largest  at  the  tip, 
just  as  for  the  ogive. 

4.3.1  Cone-sphere  Bistatic  RCS  Results  for  0.869  GHz,  Sphere-Cap  Incidence 

Two  bistatic  tests  at  0.869  GHz  are  discussed  for  the  cone-sphere.  For  the  first  group  of 
tests,  CSIX,  the  sinusoid  incident  wave  is  incident  at  0°,  or  sphere-cap  incidence.  The  second 
group  of  tests,  CS2X,  is  for  an  incident  wave  tip-on.  Two  different  grids  were  generated  for  the 
cone-sphere.  The  grid  size,  (50-73-45),  is  identical  for  both  grids,  but  the  grid  spacing  in  the 
radial  direction  is  different  at  the  surface.  For  the  first  grid,  the  spacing  is  200  cells/>.  at  the 
surface  and  for  the  second  grid,  the  spacing  is  400  cells/X,  at  the  surface. 

For  the  group  of  tests  with  sphere-cap  incidence,  both  grids  were  tested  along  with  two 
different  grid  sizes.  In  Figure  40,  the  results  using  a  coarse  grid  with  200  cells/?i  spacing  is 
shown.  Figure  40  shows  that  the  forward  scattering  from  the  tip  is  not  accurate.  This  proves  that 
this  grid  spacing  at  this  frequency  is  not  sufficient.  The  diffraction  from  the  tip  is  not  correctly 
computed  with  this  grid  spacing.  A  finer  grid  in  the  theta  and  phi  directions  also  does  not 
improve  the  results  in  the  forward  scattering  direction  as  shown  with  test  CSla.  An  increase  in 
the  spacing  in  the  radial  direction,  tests  CSlc  and  CSld,  as  shown  in  Figure  41  greatly  improves 
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RCS  (dBsm) 


■FVTD1-CS1  a  (50-73-45) 
FVTD2-CS1b  (50-39-35) 
■Moment  Method  (CICERO) 


20  40  60  80  100  120  140  160  1£ 


Theta  (Degrees) 


Figure  40:  Cone-Sphere  Bistatic  RCS,  0.869  GHz,  W,  Sphere-Cap  Incidence,  Coarse  Grid 


. FVTD1-CS1 0(50-73-45) 

- FVTD2-CS1d  (50-39-35) 

- Moment  Method  (CICERO) 


0  20  40  60  80  100  120  140  160  1 

Theta  (Degrees) 

Figure  41:  Cone-Sphere  Bistatic  RCS,  0.869  GHz,  W,  Sphere-Cap  Incidence,  Fine  Grid 


the  forward  scattering  region.  The  surface  GPD  can  be  decreased  to  22.0  cells/A.,  but  the  GPD  at 
the  surface  in  the  radial  direction  needs  to  be  high  at  this  frequency  (400  cells/X,).  The  results  for 
test  CSlc  are  slightly  closer  (<0.5  dB)  to  the  MoM  results  than  the  CSld  results.  The  small 
difference  can  be  considered  negligible;  therefore,  grid  convergence  can  be  assumed  for  the 
smaller  grid  due  to  the  computer  run  time  saved.  The  computer  run  time  saved  with  the  smaller 
grid  is  approximately  a  factor  of  seven  (22900  seconds  (CSlc)  vs.  3010  seconds  (CSld)).  Due  to 
the  time  step  increase  and  the  lower  mn  time,  test  CSld  is  considered  sufficient  and  provides  a 
very  close  computation  to  that  obtained  with  the  finer  grid  and  to  the  MoM  RCS  results. 

In  Figure  42,  the  difference  between  test  CSlc  and  the  MoM  RCS  results  are  shown  after 
the  phase  is  taken  into  account  as  discussed  in  Chapter  3.  A  difference  of  approximately  0.1  dB 
is  shown  for  the  first  40  degrees  in  the  backscatter  direction.  This  difference  is  considered 
negligible  in  the  calculations  of  RCS  data.  The  phase  difference  for  this  bistatic  data  is  not  as 
critical  as  for  the  ogive  RCS  data.  No  deep  nulls,  which  can  provide  erroneous  error 
calculations,  occur  in  the  results.  As  can  be  seen  in  the  plot,  no  more  than  0.98  dB  separates  the 
FVTD  RCS  results  from  the  MoM  RCS  results.  The  greatest  difference  occurs  in  the  forward 
scatter  region  from  the  tip.  Again,  the  tip  is  the  most  difficult  geometric  characteristic  of  the 
cone-sphere  to  accurately  compute  RCS  data.  For  the  ogive,  the  most  significant  errors  occurred 
at  the  tips.  The  same  trend  is  seen  here  with  the  cone-sphere. 

The  cross-correlation  between  the  FVTD  (CSlc)  and  MoM  RCS  results  is  shown  in 
Figure  43.  The  maximum  correlation  occurs  at  180°  which  illustrates  the  excellent  correlation 
between  the  two  sets  of  data.  The  steep  slope  of  the  cross-correlation  curve  signifies  that  no 
phase  shift  exists.  Phase  shifts  tend  to  level  the  curve  at  the  peak  or  shift  the  peak  to  either  side 
of  180°.  The  Fourier  transforms  of  both  RCS  curves  illustrate  the  excellent  agreement  between 
the  FVTD  and  MoM  RCS  results. 
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Figure  43:  Cone-Sphere  Bistatic  RCS,  0.869  GHz,  W,  Cross-Correlation  ofRCS  Data 


Figure  44:  Cone-Sphere  Bistatic  RCS,  0.869  GHz,  HH,  Sphere-Cap  Incidence,  Fine  Grid 


Similar  results  were  obtained  for  the  cone-sphere  for  0.869  GHz  for  HH  polarization. 
Figure  44  plots  the  FVTD  results  for  HH  against  MoM  RCS  data.  The  greatest  difference  of  2.5 
dB  occurs  in  the  null  at  143°.  The  FVTD  results  for  HH  miss  the  nulls  by  several  dB  for  some 
bistatic  plots.  The  same  result  was  seen  for  the  ogive  at  1.18  GHz  for  HH  in  Figure  19. 

4.3.2  Cone-Sphere  Bistatic  RCS  Results  for  0.869  GHz,  Tip-On  Incidence 

FVTD  RCS  data  for  tip-on  incidence  were  obtained  in  test  OG2X.  The  incidence  wave 
propagates  toward  the  cone-sphere  from  180°.  The  accuracy  of  the  bistatic  results  is  dependent 
on  the  correct  calculation  of  creeping  and  traveling  waves  and  diffraction  at  the  tip.  The 
creeping  waves  travel  behind  the  sphere  and  the  creeping  waves  travel  along  the  smooth  surface 
of  the  cone.  The  VV  and  HH  polarization  cases  are  shown  in  Figures  45  and  46,  respectively. 
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Figure  45:  Cone-Sphere  Bistatic  RCS,  0.869  GHz,  W,  Tip-On  Incidence 
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The  finer  grid,  CS2b,  produces  much  better  results  than  the  grid  with  the  coarser  spacing 
(CS2a)  at  the  surface.  The  forward  scatter  results  (0=0°)  are  only  accurate  if  the  fine  grid  is  used 
(CS2b).  Recall  that  the  fine  grid  has  a  spacing  of  400  cells/X  at  the  surface.  The  forward  scatter 
region  is  dependent  on  the  accurate  calculation  of  creeping  waves  as  shown  in  the  plot. 

4.3.3  Cone-sphere  Monostatic  RCS  Results  for  0.869  GHz 

The  monostatic  RCS  data  for  the  cone-sphere,  HH  polarization,  is  plotted  in  Figure  47. 
The  FVTD  and  MoM  RCS  results  match  each  other  much  closer  than  they  match  the  empirical 
data.  The  author  of  reference  [55]  states  that  errors  exist  in  the  experimental  data,  especially  in 
the  forward  scatter  area  (120°- 180°).  The  MoM  and  FVTD  results  are  almost  identical  for  every 
location  except  for  several  of  the  bistatic-to-monostatic  approximation  junctions.  The  agreement 
between  the  techniques  suggest  that  the  empirical  data  is  not  correct  from  120°- 180°. 


Figure  47:  Cone-Sphere  Monostatic  RCS,  0.869  GHz,  HH 
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Figure  48:  Cone-Sphere  Monostatic  RCS,  0.869  GHz,  W 

Monostatic  RCS  results  for  VV  polarization  for  the  cone-sphere  are  shown  in  Figure  48. 
The  I^TD  results  closely  match  the  experimental  and  MoM  RCS  results.  For  certain  regions, 
the  FVTD  results  are  closer  to  the  MoM  results  than  the  experimental  results.  The  small 
fluctuation  in  the  FVTD  curve  at  135°  shows  the  limitation  of  the  bistatic-to-monostatic 
approximation.  As  the  monostatic  curve  obtains  more  oscillations,  as  for  higher  frequencies,  the 
monostatic  simulations  every  10°  may  not  be  sufficient.  The  simulations  may  have  to  be 
performed  every  5°  to  accurately  obtain  the  oscillations  in  the  curve. 

The  FVTD  results  for  the  cone-sphere  at  0.869  GHz  using  a  sinusoid  incident  wave  were 
excellent  when  compared  to  MoM  results.  The  results  illustrate  that  for  electrically  small 
objects,  the  FVTD  algorithm  correctly  predicts  diffraction,  creeping  waves,  and  traveling  waves. 
The  following  results  are  for  a  higher  frequency,  3.0  GHz. 
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4.3.4  Cone-sphere  Bistatic  RCS  Results  for  3.0  GHz 

In  this  section,  the  bistatic  RCS  results  for  the  cone-sphere  for  3.0  GHz  are  presented. 
Test  CS4a  and  CS5a  are  the  sphere-cap  and  tip-on  incidence  tests,  respectively.  The  length  of 
the  cone-sphere  is  6.9  wavelengths  at  this  specific  frequency.  The  cone-sphere  at  this  frequency 
is  in  the  optical  region,  just  as  the  ogive  was  for  9.0  GHz.  The  electromagnetic  phenomena,  such 
as  diffraction,  creeping  waves,  and  traveling  waves,  are  local  and  the  test  time,  in  periods,  is  not 
as  long  as  for  lower  frequencies. 


Figure  49:  Cone-Sphere  Bistatic  RCS,  3.0  GHz,  HH,  Sphere-Cap  Incidence 


The  HH  polarization  case  for  test  CS4a  is  shown  in  Figure  49.  The  FVTD  RCS  results 
are  plotted  against  MoM  CICERO  data.  The  accuracy  of  the  VV  polarization  results  was  almost 
identical.  The  MoM  and  FVTD  RCS  data  differ  by  no  more  than  1.4  dB.  The  largest  difference 
occurs  in  the  forward  scatter  direction  at  the  tip  (0=180°)  and  in  the  nulls. 
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Figure  50:  Cone-Sphere  Bistatic  RCS,  3.0  GHz,  W,  Tip-On  Incidence 

Figure  50  is  the  W  polarization  result  for  tip-on  incidence  (CS5a).  The  FVTD  results 
are  almost  identical  to  the  MoM  RCS  results.  Small  differences  of  less  than  1.2  dB  occur  in 
several  of  the  oscillations  but  are  negligible.  The  FVTD  RCS  data  for  the  forward  scatter  region, 
from  the  sphere-cap  (9=0°),  differ  by  1 .0  dB  from  the  MoM  data.  The  FVTD  data  for  the 
backscatter  region  (0=180°)  from  the  tip  differ  by  1.2  dB  from  the  MoM  RCS  results.  As  seen 
with  the  ogive,  errors  in  the  FVTD  RCS  calculations  first  occur  at  the  tips.  An  increase  in  the 
GPD  will  better  the  RCS  results  but  the  computer  run  time  will  be  longer.  The  amount  of  error 
permissible  determines  the  shortest  run  time  and  GPD  required  for  the  tip  of  a  scattering  object. 

To  complete  the  analysis  of  the  cone-sphere,  the  same  metrics  are  calculated  for  the 
cone-sphere  as  for  the  ogive  to  illustrate  the  excellent  RCS  results  obtained  using  FVTD.  The 
largest  difference,  in  dB,  and  the  mean-square  error  (MSE)  are  metrics  for  presenting  the 
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difference  and  error  between  FVTD  and  MoM  RCS  results.  Table  8  lists  the  largest  dB 
difference  and  MSB  for  each  polarization  of  each  cone-sphere  test. 


Table  8:  MSE  and  Comparisons  for  the  Cone-Sphere  Tests 


Test  Number 

Difference, 
HH  (dB) 

MSE,HH 

Difference, 
VV  (dB) 

MSE,  VV 

CSla 

3.0 

3.11E-5 

2.3 

2.36E-5 

CSlb 

3.2 

3.77E-5 

2.4 

2.76E-5 

CSlc 

1.2 

2.80E-6 

1.0 

2.16E-6 

CSld 

1.6 

5.54E-6 

1.2 

4.14E-6 

CS2a 

2.4 

2.99E-5 

2.0 

3.09E-5 

CS2b 

1.3 

5.27E-5 

1.0 

4.53E-6 

CS3X 

0.5 

N/A 

N/A 

N/A 

CS4a 

1.4 

3.54E-4 

N/A 

N/A 

CS5a 

N/A 

N/A 

2.1 

2.40E-4 

The  largest  difference,  in  dB,  between  bistatic  FYTD  and  MoM  results  for  tests  CSIX- 
CS2X  is  3.2  but  this  was  for  the  grid  with  the  coarse  spacing  (200  cells/X.)  at  the  surface.  For  the 
results  using  the  finer  grid  at  the  surface  (400  cells/X.),  the  greatest  error  is  1.6  dB.  For  the 
monostatic  results,  the  largest  difference  for  HH  is  0.5  dB.  The  greatest  error  for  the  3.0  GHz 
tests  is  2.1  dB.  The  errors  for  VV  are  less  than  for  HH  just  as  for  the  ogive  RCS  results.  Test 
CSlb  had  the  largest  MSE  and  test  CSlc  had  the  smallest  MSE.  Based  on  the  plots,  test  CS4a 
visually  had  the  smallest  error  and  the  MSE  data  confirm  this  observation.  The  metrics  listed  in 
the  table  illustrate  the  excellent  FVTD  results  obtained  for  the  cone-sphere. 

4.4  Summary 

The  electromagnetic  scattering  and  RCS  results  for  the  EMCC-defined  RCS  test  bodies, 
ogive  and  cone-sphere,  were  presented.  Bistatic  and  monostatic  RCS  results  were  compared  to 


87 


MoM  and  empirical  RCS  results.  Grid  convergent  tests  were  completed  which  analyzed  the 
variation  of  the  size  of  the  grid  in  each  coordinate  direction  (R,0,<|)).  Grid  point  densities  (GPD) 
were  studied  to  arrive  at  an  optimum  GPD  in  each  direction. 

For  the  ogive,  the  bistatic  results  at  1.18  GHz  required  a  grid  spacing  of  510  cells/A<  at 
the  surface  in  the  radial  direction.  The  surface  spacing  in  the  theta  and  phi  directions  had  to  be  at 
least  22  cells/X.  The  bistatic-to-monostatic  approximation  gave  monostatic  results  for  the  ogive 
that  differed  from  MoM  RCS  data  by  no  more  than  2.5  dB  for  HH  polarization.  The  FVTD  RCS 
results  for  VV  were  similar  when  compared  to  MoM  and  experimental  RCS  data. 

The  ogive  bistatic  RCS  for  9.0  GHz  differed  by  no  more  than  2.4  dB.  The  largest 
differences  occurred  in  the  forward  and  backscatter  regions.  The  electrically  larger  ogive  at  this 
frequency  required  a  lower  PEC  surface  boundary  condition  GPD  of  200  cellsA  and  a  surface 
GPD  of  15.2-18.8  cells/^.  The  simulation  time  was  longer  but  the  number  of  periods  required 
for  convergence  was  less  than  for  lower  frequencies. 

A  Gaussian  pulse  and  sinusoid  incident  wave  produced  almost  identical  results  for  1 .0, 
2.0,  3.0,  and  4.0  GHz.  Computational  issues  such  as  numerical  dispersion  and  the  distance  of  the 
outer  boundary  from  the  scatterer  for  different  frequencies  were  critical  when  using  the  Gaussian 
pulse  incident  wave.  The  Gaussian  pulse  incident  wave  provides  results  for  multiple  frequencies 
but  the  simulation  time  is  longer  than  for  a  test  using  a  sinusoid  incident  wave. 

The  bistatic  and  monostatic  cone-sphere  results  confirmed  the  accuracy  and  grid 
requirements  for  the  ogive.  For  0.869  GHz,  the  surface  grid  spacing  required  was  22-25  cells/X. 
The  spacing  in  the  radial  direction  had  to  be  400  cells/X  to  accurately  compute  diffraction  and 
traveling  and  creeping  waves.  The  results  differed  by  no  more  than  1.6  dB  from  the  MoM  results 
and  0.5  dB  from  the  empirical  results.  The  bistatic  RCS  for  3.0  GHz  differed  by  no  more  than 
2.1  dB  from  the  MoM  results. 
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5  Conclusions  and  Recommendations 


5.1  Conclusions 

The  objectives  of  this  research  were  to 

•  modify  the  Wright  Lab  FVTD  code  to  analyze  the  electromagnetic  scattering  from 
PEC  three-dimensional  objects 

•  validate  the  characteristic-based  FVTD  formulation  by  using  the  modified  code  to 
analyze  the  electromagnetic  scattering  from  the  three-dimensional  objects,  ogive  and 
cone- sphere,  and  compare  the  FVTD  results  to  MoM  results  and  empirical  data 
published  by  the  EMCC. 

Both  objectives  were  met.  The  original  code  was  modified  to  analyze  scattering  from  closed- 
surface  perfect  electric  conductor  (PEC)  3-D  objects  using  either  a  Gaussian  pulse  or  a  sinusoid 
incident  wave.  The  specification  of  the  direction  and  polarization  of  the  incident  wave  was 
added  to  give  monostatic  and  bistatic  results.  An  RCS  convergence  check  was  also  programmed, 
used  with  a  sinusoid,  which  ends  the  simulation  after  the  transients  diminish  and  the  bistatic  RCS 
values  are  within  0. 1  dB  of  the  RCS  values  calculated  during  the  previous  period.  A  threshold 
check  used  with  a  Gaussian  pulse  was  programmed  to  end  the  simulation  once  the  scattered  field 
is  140  dB  below  the  peak  of  the  Gaussian  pulse.  A  bistatic-to-monostatic  RCS  approximation 
saves  computer  simulation  time  by  a  factor  of  nearly  40  for  computing  monostatic  data. 

Using  the  FVTD  algorithm,  the  bistatic  and  monostatic  RCS  were  calculated  and 
compared  to  MoM  and  empirical  results  in  this  research.  A  meticulous  grid  study  was  performed 
for  the  ogive  to  determine  the  grid  parameters  required  to  obtain  accurate  results.  The 
conclusions  drawn  from  the  ogive  grid  were  applied  to  the  generation  of  the  grid  for  the  cone- 
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sphere  to  obtain  accurate  results.  The  modified  code  can  be  used  to  analyze  the  scattering  from 
other  closed-surface,  single- zone,  3-D  objects  if  an  appropriate  grid  is  generated. 

5.1.1  Ogive  RCS  Results 

The  FVTD  results  for  the  ogive  were  excellent  when  compared  to  MoM  results  and 
empirical  data.  Several  bistatic  and  monostatic  tests  were  completed  for  the  ogive  at  1.18  GHz. 
The  bistatic  tests  showed  that  a  grid  point  density  (GPD)  on  the  surface  of  approximately  22-32 
cells/X  produced  the  best  results.  Shankar  reports  a  GPD  requirement  of  30-50  cells/X,  for  objects 
with  edges  or  tips,  like  the  ogive,  for  his  second-order  accurate  algorithm  [42].  The  lower  GPD 
requirement  for  Shang’s  Wright  Lab  fourth-order  accurate  FVTD  code  is  consistent  with  the 
order  of  accuracy  of  the  algorithms.  It  was  critical  for  the  spacing  in  the  radial  direction  to  be 
close  to  400-510  cells/A,  to  obtain  accurate  results  for  lower  frequencies  but  could  be  reduced  to 
200  cells/X  for  higher  frequencies  (9.0  GHz).  The  bistatic  tests  for  1.18  GHz  differed  from  the 
MoM  results  by  no  more  than  3.0  dB  for  HH  polarization  and  1.1  dB  for  VV  polarization.  The 
FVTD  calculations  for  the  monostatic  tests  were  compared  to  empirical  results  in  addition  to 
MoM  results.  The  FVTD  results  are  within  2.5  dB  of  the  MoM  monostatic  values  and  within  3.1 
dB  of  the  empirical  results. 

The  ogive  is  one  wavelength  long  at  1.18  GHz.  The  small  electrical  size  of  the  ogive  at 
this  frequency  dictates  that  electromagnetic  phenomenon  such  as  traveling  waves  have  to  be 
calculated  accurately  (<3.0  dB).  Diffraction  from  the  tip  of  the  ogive  also  contributes  to  the 
RCS.  The  small  spacing  in  the  radial  direction  (approximately  400-510  cells/)i.)  is  required  to 
accurately  compute  the  electromagnetic  scattering  due  to  these  interactions.  As  the  GPD 
decreases,  the  errors  in  the  RCS  occur  first  in  the  backscatter  and  forward  scatter  direction  as 
would  be  expected  because  of  the  diffraction  at  the  tips  of  the  ogive. 
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Two  additional  tests  were  completed  for  the  ogive  at  1.18  GHz.  The  Gaussian  pulse 
incident  wave  used  for  the  two  tests  placed  the  frequency  of  interest  near  the  center  of  the  pulse 
BW  and  near  the  edge  of  the  usable  BW  of  the  pulse.  Both  tests  showed  that  the  entire  usable 
BW  of  the  pulse  will  provide  results  almost  identical  to  those  of  the  sinusoid  incident  wave. 

The  RCS  was  calculated  for  the  ogive  at  9.0  GHz.  For  the  9.0  GHz  test,  the  surface  grid 
spacing  was  decreased  to  15-19  cells/A..  The  results  were  excellent  and  the  forward  and 
backscatter  RCS  values  differed  from  MoM  RCS  results  by  no  more  than  2.4  dB.  The  grid 
spacing  in  the  radial  direction  was  200  cells/X..  As  the  electrical  size  of  the  object  increases, 
traveling  waves  and  diffraction  contribute  less  to  the  RCS.  These  phenomena  become  local  and 
the  spacing  does  not  have  to  be  as  large  to  accurately  compute  the  propagation  of  the  wave. 

The  relatively  large  GPD  required  at  the  surface  in  the  radial  direction  as  compared  to 
the  theta  and  phi  directions  is  due  to  the  geometry  in  addition  to  the  calculation  of  specific 
electromagnetic  phenomena.  The  curves  and  the  tips  of  the  ogive  are  defined  by  the  surfaces  of 
the  finite-volume  cells.  For  the  FVTD  formulation,  the  tangential  scattered  electric  field  based 
on  the  incident  field  is  implemented  at  the  cell  centers  and  not  on  the  cell  faces  or  cell  vertices. 
The  tip  of  the  ogive  appears  to  have  a  small  blunt  tip  if  only  the  cell  centers  are  viewed.  This 
characteristic  of  the  grid  and  FVTD  formulation  requires  the  cell  size  at  the  tips  of  the  ogive  to 
be  extremely  small  to  accurately  represent  the  geometry.  The  geometry  must  be  accurately 
represented  to  correctly  calculate  the  diffraction  occurring  at  the  tips. 

The  tests  for  the  ogive  at  several  different  frequencies  show  the  length  of  time  required 
to  obtain  accurate  RCS  data  for  an  object  using  a  sinusoid  incident  wave.  Sufficient  time, 
approximately  two-four  periods  times  the  electrical  length  of  the  object,  is  required  for  the 
transients  to  diminish  due  to  the  tum-on  of  the  incident  wave.  At  lower  frequencies,  such  as  1.18 
GHz,  convergence  was  not  reached  until  four  periods.  For  a  higher  frequency,  9.0  GHz,  accurate 
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data  was  obtained  for  a  length  of  time,  in  periods,  equal  to  approximately  three  times  the 
electrical  length  of  the  object. 

The  RCS  of  the  ogive  was  obtained  using  a  sinusoid  and  a  Gaussian  pulse  incident  wave. 
The  Gaussian  pulse  is  advantageous  if  multiple  frequencies  are  desired  and  the  simulation  time  is 
less  than  the  total  time  for  the  individual  sinusoid  tests.  Excellent  results  were  obtained  using  the 
Gaussian  pulse  for  4.0  GHz  and  below.  Simulation  times  and  grid  sizes  limited  the  largest 
frequency  which  could  be  studied  using  a  Gaussian  pulse.  The  range  of  frequencies  is  limited  for 
the  Gaussian  pulse  due  to  the  location  of  the  outer  edge  of  the  grid  for  each  frequency.  Errors 
were  identified  for  the  1 .0  and  2.0  GHz  in  the  forward  scatter  region  due  to  the  location  of  the 
outer  edge  of  the  grid  (0.75^  and  1.5A,  for  1.0  GHz  and  2.0  GHz,  respectively).  The  outer  edge 
must  be  2X-3X  to  minimize  the  numerical  reflections  from  the  outer  edge  of  the  grid. 

The  metrics  established  to  compare  FVTD  RCS  results  to  MoM  and  empirical  RCS  data 
revealed  excellent  correlation.  The  greatest  difference,  in  dB,  between  the  FVTD  and  MoM  RCS 
data  was  3.0  dB  and  between  empirical  and  FVTD  data,  3.1  dB.  The  errors  for  VV  were  less 
than  for  HH  due  to  the  additional  nulls  in  the  HH  curves.  High  correlation  between  FVTD  and 
MoM  RCS  results  were  computed,  but  any  phase  shifts  showed  up  as  a  shift  of  one  or  two 
degrees  in  the  cross-correlation  curve.  The  low  MSEs  (i.e.  5.90E-1 1  to  1.08E-5)  also  depicted 
the  excellent  results. 

5.1.2  Cone-Sphere  RCS  Results 

The  cone-sphere  results  helped  to  confirm  the  ogive  results.  The  cone-sphere,  two 
wavelengths  long  at  0.869  GHz,  is  electrically  small.  Accurate  results  were  obtained  for  a 
surface  grid  spacing  of  22-26  cells/A,,  but  the  radial  spacing  for  the  PEC  surface  boundary 
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condition  must  be  400  cells/?i  to  accurately  consider  diffraction,  traveling  waves,  and  creeping 
waves. 

The  bistatic  FVTD  results  do  not  differ  by  more  than  1.6  dB  from  the  MoM  results.  The 
largest  errors  occurred  at  the  tip  of  the  cone-sphere.  The  results  could  be  improved  slightly  with 
a  finer  grid,  but  the  time  step  would  be  decreased  resulting  in  a  longer  simulation  time.  The 
slight  error  is  acceptable  and  justifiable  if  the  simulation  time  is  decreased.  Excellent  monostatic 
results  were  obtained  using  FVTD  and  differ  by  no  more  than  0.5  dB  from  the  MoM  results. 

The  bistatic  RCS  data  was  also  obtained  for  the  cone-sphere  at  3.0  GHz  using  FVTD. 
Accurate  data  was  obtained  after  9  or  10  periods  for  the  6.9  wavelength  object.  The  FVTD  RCS 
calculations  were  within  2. 1  dB  of  the  MoM  RCS  data.  Accurate  results  required  a  surface  grid 
spacing  of  14-26  cells/^  and  a  radial  spacing  of  200  cells/X,  to  accurately  consider  diffraction, 
traveling  waves,  and  creeping  waves.  These  grid  point  density  requirements  confirm  the  ogive 
conclusion  that  a  lower  grid  point  density  is  needed  for  electrically  larger  objects  since 
diffraction,  traveling  waves,  and  creeping  waves  contribute  less  to  the  RCS. 

5.1.3  FVTD  Computational  Issues 

The  experience  gained  using  FVTD  to  solve  electromagnetic  scattering  problems 
revealed  that  several  important  computational  issues  need  to  be  considered  to  obtain  accurate 
data.  The  issues  arise  because  of  the  FVTD  technique  and  numerical  algorithm,  boundary 
conditions,  type  of  incident  wave,  and  the  electromagnetic  scattering  phenomena  occurring  at  the 
surface  of  the  object. 

The  numerical  algorithm  used  for  the  flux  evaluation,  van  Leer’s  kappa  scheme  (MUSCL 
approach),  is  third-order  accurate.  The  Runge-Kutta  technique  used  for  the  time  integration  is 
fourth-order  accurate.  These  techniques  give  the  code  a  potential  high  order  of  accuracy  but  the 
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accuracy  is  degraded  if  the  grid  is  poorly  constructed.  A  grid  which  is  highly  distorted  or 
incongruent  will  produce  computational  errors.  The  grids  incorporated  grid  stretching,  such  as 
hyperbolic  tangent  stretching,  but  extreme  stretching  will  degrade  the  flux  calculations.  Further 
work  could  be  done  to  calculate  the  minimum  number  of  cells  required  in  the  radial  direction. 
For  this  research,  50-83  cells  were  used  in  the  radial  direction  but  as  few  as  15-20  cells  may  be 
sufficient  if  the  grid  is  properly  constructed.  This  improvement  would  improve  simulation  times 
by  a  factor  of  two  or  three. 

The  grid  must  be  properly  constructed  to  obtain  the  best  results  from  the  first-order 
surface  boundary  condition.  The  radial  lines  at  the  surface  of  the  scatterer  should  be 
approximately  perpendicular  to  the  surface.  Several  attempts  were  made  to  construct  grids 
analytically  but  the  grid  lines  were  not  perpendicular  to  the  surface.  These  attempts  produced 
poor  results. 

Another  computational  issue  which  is  critical  for  the  completion  of  the  computer  tests  is 
the  time  step.  The  stability  of  the  system  of  differential  equations  is  related  to  the  smallest  cell 
size.  The  ogive  and  cone-sphere  geometries  produce  converging  grid  lines  at  the  tips  resulting  in 
extremely  small  cells.  This  characteristic  for  these  geometries  revealed  a  limitation  of  the 
explicit  FVTD  algorithm  and  available  computing  resources.  The  small  time  step,  the  large  grid 
size,  and  the  number  of  periods  required  to  obtain  accurate  frequency  data  results  in  a  longer 
simulation  time.  Accuracy  in  the  backscatter  and  forward  scattering  regions  must  be  sacrificed 
for  a  faster  simulation  time  and  smaller  grid.  If  a  large  number  of  periods  required  for  the  test  is 
large,  such  as  20-25  periods,  the  resulting  computer  simulation  time  is  large  on  the  Cray  (i.e.  25- 
35  hours). 

If  the  first-order  boundary  condition  was  improved,  the  size  of  the  smallest  grid  cell 
could  be  increased  and  the  resulting  simulation  time  would  decrease.  The  time  step  is  dependent 
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on  the  cell  size.  Any  research  efforts  that  will  permit  larger  cells  relative  to  the  wavelength  will 
greatly  increase  the  flexibility  of  the  explicit  algorithm  and  potential  applications. 

The  initial  transients  produced  from  the  sinusoid  incident  wave  require  that  frequency 
data  cannot  be  taken  for  several  periods  after  the  beginning  of  the  test.  The  length  of  time  is 
dependent  on  the  electrical  size  of  the  object.  The  time  required  was  found  to  be,  in  periods,  two 
to  four  times  the  electrical  length  of  the  object.  As  the  object  increases  in  electrical  size,  the 
electromagnetic  phenomena  becomes  more  localized  and  less  time,  in  periods,  is  required. 

The  electromagnetic  phenomena  which  occurs  from  the  surfaces  of  an  ogive  and  cone- 
sphere  provide  validation  for  many  computational  codes.  The  smooth  curved  surfaces,  tips,  and 
diffraction  points  are  surface  characteristics  which  can  pose  difficulties  for  computing  accurate 
scattering  results.  Based  on  the  FVTD  results  for  the  ogive  and  cone-sphere,  the  electrical  size 
of  the  object  is  critical  when  determining  grid  size  and  spacing.  For  a  small  object,  1-2X,  the  grid 
spacing  can  be  22-25  cells/X  on  the  surface.  The  spacing  in  the  radial  direction  must  be  400-500 
cells/?i.  For  an  object  which  is  electrically  larger  (7X-SX),  the  surface  grid  spacing  may  be 
reduced  to  15-19  cells/?i,  and  the  radial  grid  spacing  can  be  decreased  to  approximately  200 
cells/A..  These  findings  are  critical  for  the  expansion  of  the  code  to  studying  electrically  larger 
objects  such  as  airfoils  and  aircraft  shaped  bodies.  A  grid  must  be  generated  which  will 
incorporate  these  features  for  a  particular  frequency. 

5.2  Suggested  Areas  for  Further  Research 

The  application  of  FVTD  to  the  area  of  CEM  is  relatively  new.  The  potential  areas  of 
future  research  are  numerous  and  several  areas  are  briefly  discussed  below. 
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5.2.1  Surface  Boundary  Condition 

The  first-order  accurate  boundary  condition  by  Shang  necessitates  a  very  dense  grid  near 
the  surface  of  the  object.  The  third-order  accurate  Van  Leer’s  kappa  scheme  used  to  calculate 
the  fluxes  through  the  faces  of  the  cells  and  the  fourth-order  accurate  Runge-Kutta  scheme  used 
for  the  time  integration  requires  a  grid  density  of  approximately  15-30  cells/X,.  The  first-order 
accurate  surface  boundary  condition  requires  a  grid  density  of  at  least  200  cells/X  for  objects 
containing  tips  such  as  the  ogive  and  cone-sphere.  An  improvement  in  the  surface  boundary 
condition  would  greatly  decrease  the  grid  density  required  at  the  surface  of  the  PEC  object. 

5.2.2  Radiation  Boundary  Condition 

For  the  purpose  of  this  thesis  research,  the  compatibility  condition  implemented  by 
Shang  at  the  outer  boundary  of  the  computational  domain  was  considered  to  be  sufficient.  The 
compatibility  condition  is  not  exact  and  the  outer  boundary  must  be  placed  at  approximately 
three  wavelengths  from  the  surface  of  the  object  to  prevent  erroneous  numerical  reflections  from 
significantly  affecting  the  scattered  field  results.  A  more  accurate  radiation  boundary  condition 
similar  to  the  perfectly  matched  layer  (PML)  boundary  condition  developed  for  FDTD  would 
greatly  reduce  the  number  of  cells  in  the  computational  space.  The  PML  developed  by  Berenger 
[4]  cannot  be  directly  implemented  in  a  characteristic-based  FVTD  formulation;  however, 
several  researchers  have  developed  generalized  PML  theories  which  could  possibly  be 
implemented  in  the  FVTD  formulation  [57]. 

5.2.3  Material  Interfaces 

For  multi-layer  dielectric  surfaces  to  be  analyzed,  the  boundary  conditions  at  the 
interface  between  the  two  materials  must  be  enforced.  This  boundary  condition  corresponds  to  a 
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jump  in  the  characteristic,  or  eigenvalue,  across  the  interface.  Further  research  needs  to  be 
completed  in  this  area. 

5.2.4  Dielectric  Materials 

Dielectric  materials  are  used  as  radar  absorbing  materials  on  wings  and  the  fuselages  of 
aircraft.  The  code  requires  further  modification  to  analyze  lossy  dielectric  materials. 

5.2.5  Frequency-Dependent  and  Time-Dependent  Materials 

As  stated  in  the  first  chapter,  a  potential  advantage  of  FVTD  is  the  ability  to  solve  the 
electromagnetic  wave  propagation  through  frequency-dependent  or  dispersive  materials  and  also 
time-dependent  materials.  A  low-frequency  or  high-frequency  CEM  code  cannot  analyze  these 
particular  types  of  materials.  Future  research  with  these  materials  would  provide  an  advantage 
for  the  FVTD  code  over  other  traditional  CEM  codes  and  techniques. 

5.2.6  Anisotropic  Materials 

Code  modifications  for  anisotropic  materials  could  also  be  considered  for  further 
research.  The  constitutive  parameters  currently  can  only  be  specified  as  scalars.  The 
constitutive  parameters  for  anisotropic  materials  are  tensors  and  the  code  requires  modification 
for  these  types  of  materials.  The  generalized  PML  theories  require  the  use  of  anisotropic 
materials  in  the  PML.  Development  in  this  area  would  be  a  step  toward  implementing  a  PML  for 
the  radiation  boundary  condition. 

RAM  structures,  such  as  honeycomb-filled  structures,  can  exhibit  anisotropic  behavior 
[42].  Development  work  and  advancement  of  the  code  for  these  materials  would  greatly  improve 
the  flexibility  and  capability  of  the  code. 
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5.2.7  Multi-zoning  for  Complicated  Objects 

Objects  coated  with  dielectric  materials  or  frequency-dependent  materials  would  require 
different  grid  structures  for  each  material.  A  zone  is  required  for  each  material  or  layer. 
Complicated  objects  may  also  require  multi-zones  for  the  accurate  generation  of  the  grid. 
Research  needs  to  be  conducted  to  expand  the  code  for  multi-zones. 

5.2.8  Hybrid  Techniques 

As  discussed  previously,  the  small  time  step  results  in  large  computation  time  due  to  the 
tips  of  the  ogive  and  cone-sphere.  Shankar  has  considered  using  high-frequency  techniques  at 
edges  and  tips  where  the  diffraction  of  the  fields  can  be  accurately  computed.  For  electrically 
large  objects,  the  diffraction  can  be  considered  a  local  phenomenon.  Hybrid  techniques  could  be 
researched  to  include  techniques  such  as  the  geometrical  theory  of  diffraction  (GTD)  or  the 
uniform  theory  of  diffraction  (UTD)  to  compute  the  diffracted  fields  at  tips  and  edges  [3].  The 
time  step  would  be  increased  and  the  resulting  computation  time  would  decrease  dramatically. 

5.2.9  Multi-discipline  Applications 

Multi-discipline  applications  could  be  studied  utilizing  the  FVTD  code  after  appropriate 
modifications.  As  discussed  previously,  the  Euler  equations  of  CFD  and  the  Maxwell  equations 
of  CEM  are  both  hyperbolic  in  nature.  The  body  of  the  code  which  calculates  the  fluxes  is 
practically  identical  for  each  application.  The  difference  in  the  code  between  the  two  disciplines 
lies  in  the  boundary  conditions,  pre-processing,  and  post-processing.  CFD  applications  require 
the  calculation  of  density  and  pressure  values  of  the  flow  field  around  the  airfoil  and  drag  and 
pitching  moment  coefficient  data  [8].  With  appropriate  additions  to  the  code,  the  aerodynamics 
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of  an  airfoil  could  be  analyzed  along  with  the  electromagnetic  scattering  from  the  airfoil.  The 
trade-offs  between  the  two  disciplines  could  be  studied  using  the  FVTD  algorithm. 

5.2.10  Code  Optimization 

The  code  is  optimized  for  a  vectorized  machine,  such  as  the  Cray  90.  Simulation  times 
could  be  decreased  if  the  code  was  optimized  for  a  parallel  machine  such  as  the  IBM  SP-2  or 
Cray  YMP.  At  higher  frequencies,  the  grid  size  required  and  the  multi-zones  for  dielectric¬ 
layered  materials  will  require  a  parallel  code  for  reasonable  simulation  times.  Research  is 
currently  being  performed  to  optimize  the  code  for  parallel  computing  machines  [6,7],  but  further 
research  has  yet  to  be  completed. 

5.2.11  Summary 

In  conclusion,  the  FYTD  code  modifications  completed  in  this  thesis  research  permitted 
the  FVTD  algorithm  and  code  to  calculate  the  RCS  of  the  ogive  and  cone-sphere.  The  FVTD 
code  results  were  excellent  as  compared  to  RCS  data  obtained  from  the  Moment  Method  (MoM) 
code,  CICERO,  and  empirical  results  published  by  the  Electromagnetic  Code  Consortium 
(EMCC).  The  comparisons  validate  the  characteristic-based  FVTD  formulation  and  code  for 
electromagnetic  scattering  problems. 
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Appendix  A:  FVTD  Formulation  and  Numerical  Algorithm 


The  finite-volume  time-domain  method  is  relatively  new  to  computational 
electromagnetics  (CEM).  As  discussed  in  the  literature  review,  Shang  and  Shankar  have 
performed  the  majority  of  the  finite- volume  time-domain  (FVTD)  effort  in  CEM.  Included  in 
this  appendix  is  a  summary  of  the  characteristic-based  FVTD  formulation  and  numerical 
algorithm  implemented  by  Shang  in  the  FORTRAN  code  used  for  this  thesis  research.  The 
FVTD  formulation  is  a  generic  finite-volume  scheme;  however,  the  numerical  technique  is  one  of 
many  schemes.  The  references  by  Shang  and  Shankar  include  numerous  other  numerical 
techniques,  used  with  the  generic  finite-volume  algorithm,  to  solve  the  Maxwell  equations  in  the 
time  domain.  This  summary  is  intended  to  aide  the  reader,  without  referencing  numerous  papers, 
in  understanding  the  FVTD  formulation  and  one  specific  numerical  technique  and  how  they  are 
implemented  to  solve  the  Maxwell  equations.  It  will  also  be  helpful  to  others  who  desire  to 
perform  subsequent  FVTD  efforts  with  application  to  CEM.  Several  equations  and  paragraphs 
from  Chapter  2  are  repeated  for  convenience  and  readability. 

A.l  Maxvrell’s  Equations 

The  time-domain  Maxwell  equations,  in  differential  form,  are  shown  below  and  will  be 
used  in  the  development  of  the  electromagnetic  FVTD  equations: 


dB 

Faraday’s  Law:  V  xE  =  - 

(A.l) 

Ampere’s  Law:  'VxH  =  ^^  +  J 

dl 

(A.2) 

Gauss’s  Electric  Law:  V  •  D  =  p 

(A.3) 

Conservation  of  Magnetic  Charge:  Y  ■  B  =  0 

(A.4) 
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where  E:  Electric  field  strength  vector  (V/m) 

D\  Electric  flux  density  vector  (C/m^) 

H\  Magnetic  field  strength  vector  (A/m) 

B\  Magnetic  flux  density  vector  (Wb/m^  or  Tesla) 

/:  Electric  current  density  vector  (A/m^) 

p:  Electric  charge  density  (C/m^) 

The  field  strength  vectors  and  the  flux  density  vectors  are  related  by  the  constitutive 
parameters.  The  constitutive  parameters,  the  electric  permittivity  and  magnetic  permeability,  are 
normally  expressed  as  tensors.  However,  if  the  material  is  linear  and  isotropic  the  constitutive 
parameters  are  scalars  and  the  constitutive  relations  become 

D  =  eE  and  B  =  pH  (A.5) 

where  e:  Electric  permittivity  (F/m) 

p:  Magnetic  permeability  (H/m) 

The  four  Maxwell  equations  are  not  independent  of  each  other.  The  two  divergence 
equations  can  be  derived  from  the  two  curl  equations  using  the  conservation  of  charge 
relationship  V  J  -  -dp/dt  assuming  that  the  material  is  linear  and  isotropic.  The  two 
divergence  equations  are  not  used  in  the  FVTD  calculations  but  can  be  used  as  a  check  on  the 
predicted  field  response  [18]. 

A.2  Maxwell’s  Equations  in  Conservation  Form 

For  use  in  FVTD,  the  two  curl  Maxwell  equations  are  cast  in  conservation  form  [37]. 

The  conservation  form  is  not  required  for  the  Maxwell  equations  but  is  required  for  the  Euler 
equations  of  fluid  dynamics.  For  the  Euler  equations,  the  conservation  form  conserves  physical 
properties  such  as  energy,  mass,  and  momentum  [8].  To  take  advantage  of  the  same 
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computational  tool,  the  Maxwell  equations  are  also  cast  in  conservation  form.  To  place  the  two 
curl  Maxwell  equations  in  conservation  form,  the  curl  operations  are  carried  out,  the  constitutive 
relations  are  implemented,  and  the  result  is  given  by  [37] 


^  dF  ^  dG  ^  dH 
dt  dx  dy  dz 


where 


(A.6) 
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Equation  (A.6)  is  a  system  of  six  linear  equations  with  six  unknowns  (U  vector).  The  equations 
are  linearly  dependent;  therefore,  a  characteristic-based  scheme  is  developed  to  uncouple  the  six 
equations.  The  system  constitutes  a  hyperbolic  system  of  partial  differential  equations  and  an 
initial  value  problem  [34].  The  hyperbolic  system  possesses  real  eigenvalues  and  independent 
eigenvectors. 


A.3  Coordinate  Transformation 

To  analyze  the  scattering  of  various  shapes  including  the  ogive  and  cone-sphere  analyzed 
in  this  thesis,  a  curvilinear  coordinate  transformation  is  required  [37].  The  transformation 
defines  a  one-to-one  relationship  between  two  sets  of  temporal  and  spatially  dependent  variables. 
The  variables  Tj,  and  ^  are  used  and  are  all  functions  of  x,  y,  and  z.  Equation  (A.6)  after  a 
coordinate  transformation  becomes  [37,  38] 


dU  dF  dG  dH 

- 1 - 1 - 1 - 

dt  d^  dq  dC 


(A.7) 
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where  U  =  — 
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-  J 


F  = 


G  = 
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V  is  the  Jacobian  of  the  coordinate  transformation  and  is  given  by  [34] 
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(A.8) 


For  example,  F  is  transformed  to  F  and  is  equal  to  [34] 
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(A.9) 


G  and  H  are  the  same  as  F  ,  except  ^  is  replaced  with  T]  and  respectively.  F,  G,  and  H  are 
the  flux  vectors  and  represent  the  electric  and  magnetic  flux  densities  at  the  faces  of  each  finite- 
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volume  cell  of  the  discretized  space  in  transformed  curvilinear  coordinates.  The  finite- volume 
formulation  of  Equation  (A.7)  and  the  calculation  of  the  flux  at  each  cell  face  will  be  discussed 
next. 


A.4  Finite- Volume  Formulation 

Equation  (A.7)  is  applied  to  every  finite-volume  cell  in  the  grid.  An  integration  is 
performed  over  each  finite-volume  cell: 


f  — 111/- 


(A.  10) 


The  divergence  theorem  is  then  applied  to  the  second  integral: 


where  n:  Unit  vector  normal  to  the  surface  (^,  T),  and  ^  for  F,  G,  and  H,  respectively) 

S:  Closed  surface  bounding  the  finite  volume  (m^) 

If  the  size  of  the  cell  is  invariant  with  time  and  a  finite-volume  discretization  is  applied 
to  the  closed  surface  integral,  then  Equation  (A.  11)  becomes  [7] 

f^v  +  ^(/t5  +  c„+H^)=-yv  (A./2) 

m=I 


where 


m: 


V: 


Magnitude  of  the  F  flux  vector  through  face  m  in  the  ^  direction 
Magnitude  of  the  G  flux  vector  through  face  m  in  the  T|  direction 
Magnitude  of  the  N  flux  vector  through  face  m  in  the  ^  direction 
Index  for  the  sides  of  the  cell 

Volume  of  cell,  equal  to  the  Jacobian  of  the  coordinate  transformation 
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Equation  (A.  12)  is  applied  to  every  volumetric  cell  in  the  grid.  The  equation  is  the 
generic  FVTD  formulation.  The  components  of  U  are  physically  located  at  the  center  of  the  cell. 

The  fluxes,  F,  G,  and  H,  are  located  on  the  surfaces  of  the  cell  and  represent  the  tangential 
components  of  the  electric  and  magnetic  fields  [8]. 

To  evaluate  the  magnitude  of  the  flux  components  through  each  face,  the  magnitude  of 
the  flux  density  in  each  coordinate  direction  is  multiplied  by  the  area  of  the  face.  The  area  of 
each  face  is  evaluated  by  taking  the  dot  product  of  the  normal  and  the  surface  area  vector, 
essentially  the  magnitude  of  the  surface  area  vector  [8]: 


F„^  =  F(ndS„-^) 

(A.13) 

G„^=G(ndS,„Ti) 

(A.14) 

H„^=^(ndS„0 

(A.15) 

A  representation  of  a  finite- volume  cell  is  shown  in  Figure  A.  1 .  The  surface  area  vectors 

required  in  Equations  (A.13-A.15)  to  calculate  the  magnitude  of  the  fluxes  through  each  face  are 

calculated  as  follows,  using  face  2  as  an  example: 

ndS2  —  T|dS2  —  ^  (Res  ^  ^72) 

(A.16) 

where  T]:  Unit  vector  normal  to  face  2 

Rea:  Vector  joining  cell  vertices  6  and  3 

R72:  Vector  joining  cell  vertices  7  and  2 

Rea  and  R72  are  found  using  the  vertices  of  face  2: 

R63  =  (x6 -Xjji  +  Cy^ -yajj +(Z6  -Zjjk 

(A.17) 

R72  =  (X7  -  X2)i  +  (y7  -  y2)j  +  -  Z2)k 

(A.18) 
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Figure  A.  1:  Finite-Volume  Cell 


The  magnitudes  of  each  flux  component  in  Equations  (A.13)-(A.15)  are  calculated  using 
a  flux-splitting  method.  The  split  fluxes  correspond  to  the  magnitude  of  the  flux  associated  with 
the  positive  and  negative  eigenvalues.  The  magnitude  of  the  flux  is  the  sum  of  the  individual 
split  fluxes.  The  details  for  the  flux-splitting  will  be  discussed  in  Section  A.6. 

Equation  (A.  12)  is  a  generic  finite- volume  formulation.  One  of  numerous  techniques  can 
be  used  to  calculate  the  fluxes  and  to  perform  the  time  integration.  Shang  uses  a  Monotone 
Upstream-Centered  Scheme  for  Conservation  Laws  (MUSCL)  to  calculate  the  fluxes  on  the 
surfaces  of  the  cells  and  a  multi-stage  Runge-Kutta  method  to  solve  it  in  the  temporal  domain 
[37].  The  methods  used  to  solve  it  are  characteristic-based;  therefore,  a  brief  review  of 
applicable  linear  algebra  will  be  given  to  provide  a  foundation  for  the  discussion  of  the 
characteristic-based  scheme. 
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A.5  Eigenvalues  and  Eigenvectors 

As  review,  a  system  of  linear  equations  can  be  solved  using  linear  algebra.  The  Maxwell 
equations  in  conservation  form  are  a  complete  set  of  six  linear  equations.  The  system  of 
equations  can  be  solved  using  eigenvalue  and  eigenvector,  or  characteristic,  analysis.  The 
eigenvalues  and  eigenvectors  of  a  linear  system  of  equations  satisfy  the  equation  [49] 

(T>-Xl)x  =  0  (A.19) 


where  D:  Coefficient  matrix 
X:  Eigenvalue 
I:  Identity  matrix 

x:  Eigenvector  corresponding  to  an  eigenvalue 

The  eigenvalues  and  eigenvectors  of  the  coefficient  matrices  in  transformed  coordinates 
is  not  required  for  the  present  FVTD  formulation  [37].  Only  the  eigenvalues  and  eigenvectors  of 
the  coefficient  matrices  of  F,  G,  and  H  in  Cartesian  coordinates  are  required  because  of  the  local 
orthogonal  coordinate  system  developed  and  discussed  later.  The  flux  vectors  F,  G,  and  H,  can 
be  represented  as  F=AU,  G=BU,  and  H=CU,  respectively,  if  the  constitutive  parameters  are 
scalars.  U  is  the  same  vector  as  used  in  Equation  (A.6).  The  coefficient  matrices  A,  B,  and  C 
are  shown  below  [29]: 
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(A.20) 
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0  0  0  0  -0 
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(A.22) 


-  0  0  0  0  0 


0  0  0  0  0  0 


The  eigenvalues  of  A,  B,  and  C  are  [29] 


*  r.-  .  1  1  1  1  ^  ^ 

A  =  Diagonal] -p=r,  . — ,  — =•,  — p=,  0,  0 

yjm 


(A.25j 


The  eigenvalues  are  equal  to  the  speed  of  the  electromagnetic  wave  in  the  medium.  In  free 
space,  the  eigenvalues  are  equal  to  the  speed  of  light.  The  signs  indicate  the  direction  of 
propagation  of  the  wave.  The  eigenvalues  are  repeated  and  contain  multiplicities,  but  linearly 
independent  eigenvectors  can  still  be  found  [34]. 

The  coefficient  matrices  A,  B,  and  C  can  be  transformed  using  similarity  matrices.  The 
transformation  uncouples  the  six  equations  in  the  system  of  linear  equations.  Equation  (A.6),  so 
that  each  can  be  solved  independently.  The  transformation  has  the  form  [49] 


A=S'‘DS 


(A.24) 


where  A:  Diagonal  matrix  containing  the  eigenvalues  of  D  along  the  diagonal 
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For  the  FVTD  formulation,  split  fluxes  are  calculated  [37].  The  flux- splitting  procedure 
splits  the  fluxes  according  to  the  direction  of  propagation  [34].  The  split  fluxes  are  calculated 
using  the  similarity  matrices  and  the  corresponding  positive  or  negative  eigenvalue.  Forward 
differencing  is  used  for  the  negative  eigenvalue  and  backward  differencing  is  used  for  the 
positive  eigenvalue  [30].  The  calculation  of  the  split  fluxes  is  facilitated  with  the  use  of  a  local 
orthogonal  system. 

A.6  Flux  Evaluation  using  a  Local  Orthogonal  Coordinate  System 

In  Cartesian  coordinates,  the  split  flux  vectors  can  be  calculated  by  [37] 


F\  =SA^S''t/\ 

i+-  i+- 

2  2 

(A.28) 

F",  =SA'S"'t/‘  , 

i+-  i+- 

(A.29) 

where  Diagonal  matrix  containing  the  positive  eigenvalues  which  correspond  to  the 

flux  flowing  in  the  positive  direction  through  the  face  with  index  i+1/2 
A':  Diagonal  matrix  containing  the  negative  eigenvalues  which  correspond  to  the 

flux  flowing  in  the  negative  direction  through  the  face  with  index  i-t-1/2 
The  calculation  of  the  flux  in  the  characteristic-based  FVTD  formulation  in  transformed 
coordinates  is  greatly  facilitated  with  the  use  of  a  local  orthogonal  system  in  the  transformed 
space  (i.e.  T),  and  ^  space).  The  local  orthogonal  system  is  established  for  each  face  of  a  cell 

and  includes  the  normal  to  the  face  and  two  unit  vectors  lying  in  the  face.  Diagonalization  of  the 
coefficient  matrices  in  the  curvilinear  coordinate  system  is  not  required  [37].  The  transformation 
between  two  orthogonal  systems  is  well  known  and  permits  easy  calculation  of  the  fluxes: 

F-  =  M"' W*  =  M“‘F-  (A.30) 
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(A.31) 


F  1  =M"'(S-'A+Sf7^  +S''A‘Sf/-  )  =  F%+F;^ 

i+—  i+i.  i+i.  in —  in — 

2  2  2  2  2 

where  M"':  Transpose  of  the  matrix  containing  the  metrics  of  the  coordinate  transformation 

M:  Matrix  containing  the  metrics  of  the  coordinate  transformation  between  the 

Cartesian  coordinate  system  and  the  local  orthogonal  coordinate  system 

M:  Matrix  containing  the  metrics  of  the  coordinate  transformation  between  the 

Cartesian  coordinate  system  and  the  curvilinear  coordinate  system 

F~ :  Positive  and  negative  split  flux  through  the  face  where  the  local  orthogonal 

system  has  been  established 

^  *4“ 

F  :  Positive  and  negative  flux  through  a  face,  expressed  in  curvilinear  coordinates 
^  , 

F  ,  :  Positive  split  flux  expressed  in  curvilinear  coordinates  through  the  (i+1/2)  face 

i*i— 

2 

F"  1 :  Negative  split  flux  expressed  in  curvilinear  coordinates  through  the  (i+1/2)  face 

i+- 

2 

F  1 :  Sum  of  the  positive  and  negative  flux  expressed  in  curvilinear  coordinates 

i+- 

2 

through  the  (i+1/2)  face 

S  :  Matrix  containing  the  eigenvectors  of  the  local  orthogonal  system,  same  as  Sa  in 

Equation  (A.26) 

The  flux  which  is  evaluated  using  Equation  (A.31)  is  used  in  Equation  (A.  13).  Note  that 
the  flux  is  evaluated  at  the  face  and  the  values  of  the  U  vector  on  the  face  must  be  calculated 
from  the  cell  centers.  The  numerical  accuracy  of  the  numerical  algorithm  is  dependent  on  the 
extrapolation  procedure  used  to  calculate  the  fluxes  from  the  values  of  the  fields  at  the  cell 
centers.  The  accuracy  is  also  dependent  on  the  metrics  of  the  cells  (i.e.  surface  area  vector)  [34]. 
The  same  flux-splitting  procedure  is  used  to  calculate  the  G  and  H  flux  vectors. 
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The  values  of  the  U  vector  on  the  cell  interfaces  are  calculated  using  a  high-order 
accurate  extrapolation  scheme  developed  by  van  Leer  [53],  The  Monotone  Upstream-Centered 
Scheme  for  Conservation  Laws  (MUSCL)  is  used  to  create  a  third-order  accurate  solution  for  the 
Maxwell’s  equations.  The  fluxes  in  Equation  (A.32)  and  the  G  and  H  fluxes  are  calculated 
based  on  piece-wise  data  from  the  centers  of  adjacent  cells  {U  vector)  from  the  previous  time 
step.  Values  of  the  U  vector  are  required  at  the  cell  face  to  calculate  the  split-flux  vectors.  The 
fields  for  one  side  of  the  cell,  such  as  i  =  1/2,  are 

1  =U,+^[(\ -K)(U,  -  U,_, )  -h (1  +K)(Ui„  -  U, )]  (A.32) 

i+—  4 

2 

-7[(l+K)(f/w  -Ui)  +  (l-K)(t/i,2  -Ui,,)]  (A.33} 

‘+2  4 

where  k  and  (j)  are  parameters  which  control  the  accuracy  of  the  approximation.  To  obtain  third- 
order  accuracy  from  van  Leer’s  kappa  scheme,  k  =  1/3  and  (j)  =  1 .  The  subscripts  correspond  to 
the  cell  center  locations  referenced  to  the  current  cell. 

A.7  Time  Integration 

For  the  time  integration  of  Equation  (A.  12),  a  multi-stage  Runge-Kutta  scheme  is  used 
[38].  A  fourth-order,  four-stage  scheme  is  calculated  as  follows: 


U®  =t/®(t,(/'’) 

(A.34) 

t/'  =  U\i  +  MI2,W  +MI2-U^) 

(A.35) 

t/^  =  t/^(t  +  At  /  2,  t/" -1- At  /  2  •  [/') 

(A.36) 

[/^  =  i7\t  +  At,U"-)-At-{/^) 

(A.37) 

=W+MI6-{lf+2U^+2U^  +  U^) 

(A.38) 

where  [/":  Value  of  the  electric  and  magnetic  flux  densities  at  the  cell  centers  at  stage  m 
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For  other  orders  of  accuracy,  refer  to  [26],  The  intermediate  steps  between  Un  and  Un+i 
have  no  physical  meaning.  The  values  of  U  are  physically  located  at  the  cell  center.  V  is  found 
by  adding  the  split  fluxes  and  dividing  by  the  cell  volume  as  in  Equation  (A.  12). 


A.8  Incident  Wave 

To  obtain  the  scattered  field  from  an  object,  an  incident  field  must  be  generated  in  the 
computational  grid.  Various  waveforms  can  be  used  to  analyze  the  scattering  from  an  object.  A 
sinusoid,  sin(©t),  can  be  used  to  analyze  the  scattering  for  a  single  frequency.  However,  an 
advantage  of  the  time-domain  analysis  is  that  multiple  frequencies  can  be  analyzed 
simultaneously  with  one  computer  simulation.  A  Gaussian  pulse  is  ideal  to  use  as  an  incident 
wave  for  multiple-frequency  analysis.  The  Gaussian  pulse  contains  multiple  frequencies  and  has 
the  form  [18] 


g(t)  =  exp 


(A.39) 


where  g(t):  Gaussian  pulse  with  an  amplitude  of  unity 
to:  Center  of  Gaussian  pulse 

T :  Period/duration  of  the  Gaussian  pulse 

The  Gaussian  pulse  and  the  parameters  for  the  incident  wave  are  described  in  detail  in 
Chapter  3.  A  trade-off  exists  between  a  sinusoid  incident  wave  and  a  Gaussian  pulse.  A 
Gaussian  pulse  will  give  results  for  multiple  frequencies,  but  the  duration  of  the  simulation  can 
be  several  times  longer  than  the  test  with  the  sinusoid  incident  wave.  The  particular  incident 
wave  can  be  selected  depending  on  the  results  desired. 
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A.9  Boundary  Conditions 


A.9.1  Scatterer  Surface  Boundary  Conditions 

For  perfect  electric  conductor  (PEC)  objects,  surface  boundary  conditions  state  that  the 
tangential  components  of  the  total  electric  field  are  equal  to  zero  and  the  normal  components  of 
the  total  magnetic  flux  density  are  zero: 

nxE  =  0  (A.40) 

n  B  =  0  (A.41) 

In  the  FVTD  formulation,  the  electric  and  magnetic  fields  are  both  calculated  on  the  surface  of 
the  scatterer.  Therefore,  extrapolated  boundary  conditions  are  needed  to  solve  for  the  unknowns. 
The  extrapolated  boundary  conditions  which  are  introduced  to  solve  for  the  total  fields  at  the 
surface  of  the  scatterer  are  given  by  [38] 

n-V(lnx(H, -H2)I)  =  0  (A.42) 

n-V(n-(D, -D2))  =  0  (A.43) 

where  n  x  is  -  electric  current  density  with  finite  jump  of  constant  value  [38] 

n  (Dj  -D2):  p  -  electric  charge  density  with  finite  jump  of  constant  value  [38] 

The  current  and  charge  densities  are  assumed  to  be  finite  jumps  and  piece-wise  continuous. 

A.9.2  Radiation  Boundary  Condition 

In  a  computational  space,  the  truncated  space  will  produce  erroneous  errors.  The 
truncated  space  produces  reflections  which  propagate  through  the  grid.  Shang  uses  a 
compatibility  condition  which  sets  the  incoming  flux  component  equal  to  zero  at  the  outer  edge 
of  the  computational  grid. 
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Ideally,  the  one-dimensional  characteristic  boundary  condition  is  exact  for  wave  motion 
in  one  coordinate  direction.  For  multiple  dimensions,  the  boundary  condition  reduces  to  first- 
order  accuracy  depending  on  the  wave  direction.  Because  of  the  low  order  of  accuracy,  the  outer 
edge  of  the  grid  must  be  placed  two  to  three  wavelengths  from  the  scatterer  surface  to  minimize 
the  reflections  which  contaminate  the  scattering  computations. 


A.IO  Fourier  Transform 

A  Fourier  Transform  (FT)  can  be  used  to  obtain  frequency  data  from  the  time-domain 
scattered  field  results.  The  FT  takes  the  following  form  [23]: 


X((0)  =  ^^x(n)e-j“'^' 

^  n=0 


{A.44) 


where  n:  Index  for  time  increments 

N :  Maximum  time  increment  for  the  FT 

co:  Radian  frequency 

x(n):  Total  field  component  calculated  at  time  n 

X(co):  Amplitude  of  the  scattered  electric  or  magnetic  field  for  radian  frequency  co 
The  FT  results  are  used  in  the  near-to-far  field  transformation.  The  transformation 
requires  sinusoidal  or  phasor  values  for  the  total  fields.  The  summation  in  Equation  (A.45)  is 
performed  for  each  frequency  at  the  end  of  each  time  step  as  an  accumulating  sum. 


A.ll  Near-to-Far  Field  Transformation 

To  obtain  far-field  parameters,  such  as  the  radar  cross  section  (RCS),  a  near-field  to  far- 
field  transformation  is  required.  The  transformation  requires  equivalent  electric  and  magnetic 
currents  and  charges  on  a  virtual  surface  surrounding  the  object  of  interest.  A  review  of  the 
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surface  equivalence  theorem  will  be  given,  and  then  the  near-field  to  far-field  transformation  will 
be  discussed  in  detail. 


A.11.1  Surface  Equivalence  Theorem 

The  surface  equivalence  theorem,  or  Huy  gen’s  principle,  places  a  virtual  surface  around 
an  object  or  sources  and  replaces  the  electric  and  magnetic  fields  on  the  boundary  of  the  virtual 
surface  with  equivalent  electric  or  magnetic  currents  and  charges  which  satisfy  the  boundary 
conditions  [3],  The  currents  are  selected  so  that  the  fields  inside  the  surface  are  zero  and  the 
fields  outside  are  equivalent  to  the  fields  produced  by  the  sources  on  the  surface  or  equivalent  to 
the  scattered  fields.  The  currents  are  calculated  from  the  tangential  fields  over  the  virtual  surface 
and  the  equivalent  charges  are  calculated  from  the  normal  fields  at  the  surface  [51]: 


Js  =  n  X  H 

(A.45) 

Ms  =  -  n  X  E 

(A.46) 

Pes=n  H 

(A.47) 

Pms  =  nE 

(A.48) 

where 


n:  Outward  normal  of  the  virtual  surface 

E,  H:  Phaser  representations  of  the  total  electric  and  magnetic  fields 
Js:  Equivalent  electric  current  density  tangential  to  virtual  surface 

Ms!  Equivalent  magnetic  current  density  tangential  to  virtual  surface 
Pes:  Equivalent  electric  charge  on  the  virtual  surface 

Pms:  Equivalent  magnetic  charge  on  the  virtual  surface 

The  equivalent  currents  and  charges  are  used  in  the  near-field  to  far-field  transformation. 
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A.11.2  Transformation 


The  equivalent  currents  and  charges  are  used  to  find  the  far  fields;  however,  the 
transformation  requires  time-harmonic  quantities.  The  equivalent  phasor  currents  and  charges 
are  calculated  from  the  scattered  fields  by  using  a  FT  and  the  equivalence  theorem  as  discussed 
in  the  previous  section  [51].  Once  the  currents  and  charges  are  found,  the  following 
transformation  computes  the  far-field  quantities  [12]  for  the  scattered  field: 


E*  = 


exp(jkR)  jk  rf 
R  471  Ws 


— (n  X  H)  -  (n  X  E)  X  r  -  (n  •  E)r 


exp(-jkr  •  R'  )dS' 


{A.49) 


H® 


exp(jkR)  jk 


R 


47ctt' 


I 

a:; 


(n  X  E)  -I-  (n  X  H)  X  r  +  (n  •  H)r 


exp(-jkr  •  R  )dS'  (A.50) 


where  E^H^  Far-field  scattered  electric  and  magnetic  fields 
k:  Wave  number 

n;  Unit  vector  normal  to  the  virtual  surface 

r:  Unit  vector  in  the  direction  of  the  far-field  observation  point 

R:  Vector  from  the  origin  to  the  observation  point 

R':  Vector  from  the  origin  to  a  point  on  the  virtual  surface 

S':  Closed  virtual  surface  over  which  the  equivalent  surface  charges  and  currents  are 

integrated 

The  primed  coordinates  refer  to  an  integration  over  the  virtual  surface  or  a  point  on  the 
virtual  surface  as  shown  in  Figure  A.2.  The  observation  point  is  in  the  far  field  and  the  various 
vectors  describing  the  location  of  the  scatterer  and  observation  point  are  shown  in  the 
illustration. 
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Observation 


Figure  A.2:  Virtual  Surface  and  Far-Field  Sketch 


A.12  RCS  Calculations 

The  RCS  can  be  obtained  for  the  backscatter  direction  or  the  bistatic  direction.  The  RCS 
is  calculated  by  using 


a 


3-D 


=  lim  47iR^ 

R — 


(A.51) 


where  E‘:  Incident  electric  field  (sinusoid  or  Gaussian  pulse) 

E®:  Scattered  electric  field 

R:  Distance  to  the  observation  point  (Magnitude  of  R  in  Equation  (A.50)) 
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The  magnitude  of  the  scattered  electric  field  is  calculated  using  Equation  (A.50).  The  RCS  can 
also  be  determined  from  the  scattered  magnetic  field  in  the  same  way  using  Equation  (A.51)  and 
the  equivalent  equation  of  (A.52)  for  the  magnetic  field. 

A.13  Summary 

The  FVTD  formulation  calculates  the  scattering  and  the  RCS  from  PEC  objects  using  a 
scattered-field  formulation  of  the  time-dependent  Maxwell  equations.  A  characteristic-based 
FVTD  formulation  is  used  to  solve  the  Maxwell  equations.  The  FVTD  formulation  uses  a  van 
Leer’s  kappa  scheme  for  the  flux  evaluation  and  a  Runge-Kutta  multi-stage  scheme  for  the  time 
integration.  The  far-field  scattering  results  such  as  the  RCS  are  obtained  from  the 
electromagnetic  fields  subsequent  to  a  Fourier  transform  and  a  near-to-far  field  transformation. 
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Appendix  B:  Finite- Volume  Time-Domain  FORTRAN  Code 


This  appendix  contains  the  code  listings  and  descriptions  of  the  FORTRAN  finite- 
volume  time-domain  (FVTD)  code  modifications.  Also  included  in  this  appendix  is  a  MATLAB 
code  for  comparing  FVTD  results  to  Moment  Method  results  and  a  MATLAB  code  for 
generating  a  movie  of  contour  slices  of  the  scattered  electric  field  from  a  cone-sphere. 

B.l  FVTD  FORTRAN  Code  Outline 

The  basic  structure  of  the  FORTRAN  77  code  is  shown  in  Figure  B.l.  The  code  closely 
follows  the  outline  in  Appendix  A  which  describes  the  procedures  required  to  solve  an 
electromagnetic  scattering  problem  using  FVTD.  Apart  from  the  time  loop  and  the  loops 
required  to  calculate  the  fluxes,  the  code  follows  a  sequential  process  to  solve  for  the  RCS  of  a  3- 
D  object. 

At  the  beginning  of  the  code,  the  input  data  and  the  grid  data  are  read.  Initial  parameters 
for  the  incident  wave  and  metrics  of  the  finite-volume  cells  are  computed.  The  time  loop  is  then 
entered  and  for  each  time  step,  the  F,  G,  and  H  flux  vectors  are  calculated  in  separate 
subroutines.  Boundary  conditions  are  implemented  at  the  beginning  of  each  time  step.  The 
fluxes  are  computed  for  the  faces  of  the  finite-volume  cells  as  described  in  Appendix  A.  The  U 
vector,  physically  located  at  the  center  of  the  cells,  is  computed  from  the  fluxes  at  the  end  of 
each  time  step.  Once  the  time  loop  is  completed,  the  Fourier  transform  and  near-to  far  field 
transformation  are  calculated,  and  the  RCS  values  are  computed  from  the  field  data. 

As  shown  in  Figure  B.l,  the  Runge-Kutta  loop  is  completed  four  times  for  every  time 
step.  The  fourth-order  accurate  code  requires  the  fluxes  (electric  and  magnetic  fields)  to  be 
computed  four  times  for  every  time  step  to  obtain  the  fourth-order  accuracy.  The  extra 
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computational  time  required  would  be  excessive  except  the  higher-order  accurate  code  permits  a 
much  larger  CFL  value  resulting  in  a  larger  time  step  as  compared  to  a  lower-order  accurate 
code. 

B.2  FVTD  Code  Listings  of  ModiHcations 

The  original  code  contained  approximately  1450  lines  of  code  and  the  modified  code 
contains  approximately  2350  lines  of  code.  The  modifications  include  the  following  subtopics; 

•  input  file  (FVTD.DAT)  for  specifying  grid  size,  geometry  type,  incident  wave  type, 
material  parameters,  stability  parameters,  and  Fourier  transform  variables 

•  input  grid  file  with  the  option  to  add  or  delete  grid  points  depending  on  the 
frequency  of  the  incident  wave  and  the  desired  accuracy 

•  incident  wave  type,  polarization,  and  direction 

•  RCS  convergence  check  and  threshold  check  used  in  conjunction  with  a  sinusoid  and 
Gaussian  pulse  incident  wave,  respectively 

•  bistatic-to-monostatic  approximation  to  obtain  monostatic  values  from  bistatic 
calculations 

The  code  modifications  are  described  and  listed  below. 

B.2.1  Input  File 

The  input  file,  FVTD.DAT,  for  the  FVTD.F  code  is  listed  below.  The  input  parameters 
are  described  in  the  file.  Parameters  for  the  grid  size,  geometry,  incident  wave,  material, 
stability,  and  Fourier  transform  can  be  specified  in  the  input  file.  Most  of  the  variables  are  self- 
explanatory  but  several  clarifications  are  required. 
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*******************  fvtD.DAT:  input  FILE  FOR  FVTD.F  ******************* 
INPUT  VARIABLES: 

GRID  VARIABLES: 

IL:  MAX  X/XI  INDEX,  CORRESPONDS  TO  NUMBER  OF  GRID  POINTS 

JL:  MAX  Y/ETA  INDEX,  CORRESPONDS  TO  NUMBER  OF  GRID  POINTS 

KL:  MAX  Z/ZETA  INDEX,  CORRESPONDS  TO  NUMBER  OF  GRID  POINTS 

GEOMETRY  VARIABLES 

IGEO:  INDICATOR  FOR  TYPE  OF  GEOMETRY 

"0":  SPHERE  (GRID  GENERATED  IN  FVTD  CODE) 

"1":  OGIVE  (GRID  FILE,  OGIVE. GRD,  REQUIRED) 

"2":  CONE-SPHERE  (GRID  FILE,  CONESPH.GRD,  REQUIRED) 

"3":  MISCELLANEOUS  GEOMETRY  (GRID  FILE  REQUIRED) 

RO:  DISTANCE  OF  OUTER  BOUNDARY  BEYOND  SCATTERER  SURFACE 

(IN  TERMS  OF  WAVELENGTHS) 

NGPD:  TIP  GRID  POINT  DENSITY  IF  REQUIRED  ( POINTS /WAVELENGTH) 

INCIDENT  WAVE  VARIABLES: 

AMP:  PEAK  AMPLITUDE  OF  THE  INCIDENT  WAVE 

PHI:  PHI  DIRECTION  OF  INCIDENT  WAVE  (SPHERICAL  COORDINATES) 

THI:  THETA  DIRECTION  OF  INCIDENT  WAVE  (SPHERICAL  COORDINATES) 

EPHI:  E  FIELD  COMPONENT  (POLARIZATION)  IN  THE  PHI  DIRECTION 
ETHI:  E  FIELD  COMPONENT  (POLARIZATION)  IN  THE  THETA  DIRECTION 
NOTE:  SQRT(EPHI**2  +  ETHI* *2)  SHOULD  BE  EQUAL  TO  ONE 

ISRC:  INDICATOR  FOR  TYPE  OF  INCIDENT  WAVE 

"0":  SINUSOID  (SPECIFY  W) 

"1":  GAUSSIAN  PULSE  (SPECIFY  GPT,  GPD,  AND  WD) 

W:  SINUSOID  RADIAN  FREQUENCY 

GPT:  GAUSSIAN  PULSE  PERIOD  (NORMALIZED:  (T*C) ) 

GPD:  GAUSSIAN  PULSE  TIME  DELAY  (CENTER  OF  PULSE,  IN  TERMS  OF  GPT) 

WD:  GAUSSIAN  PULSE  RADIAN  FREQUENCY  DELTA  (FREQUENCY  INCREMENT) 

MATERIAL  VARIABLES: 

REP:  RELATIVE  EPSILON,  ELECTRIC  PERMITTIVITY 

RMU:  RELATIVE  MU,  MAGNETIC  PERMEABILITY 

STABILITY  VARIABLES: 

CFL:  COURANT-FRIEDRICHS-LEWY  NUMBER,  CONTROLS  STABILITY 

FOURIER  TRANSFORM  VARIABLES: 

ICON:  CONVERGENCE  INDICATOR  FOR  SINUSOID  INCIDENT  WAVE 

" 0 " ;  DO  NOT  USE  CONVERGENCE  CHECK 

"1":  USE  CONVERGENCE  CHECK  (NBS,NES,AND  NEND  ARE  NOT  USED) 

NBS:  TIME  BEGINNING  POINT  FOR  THE  FT  (NUMBER  OF  PERIODS) 

NES:  TIME  ENDPOINT  FOR  THE  FT  (NUMBER  OF  PERIODS) 

NEND:  TIME  ENDPOINT  FOR  THE  FVTD  CODE  (NUMBER  OF  PERIODS) 

•k'k'k’k’k'k'k'k'k'k'k'k'k'k'kif'k'k'k'k'h'k'k'kif'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'ie’k'k'k'k'k'kic'kir'k'ic'k'k'k'k'k'k'k'k'k'kic'k-kir'k'kicir'k 


INPUT 

# 

PARAMETERS 

IL 

71 

JL 

74 

KL 

45 

IGEO 

1 

RO 

3.0 

NGPD 

30 

AMP 

1.0 

PHI 

0.0 

THI 

0.0 

EPHI 

0.0 

ETHI 

1.0 

ISRC 

0 

W 

24.17386 

GPT 

0.0182556 

GPD 

4.0 

WD 

10.4719755 
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REP  1 . 0 
RMU  1 . 0 
CFL  1 . 5 
ICON  0 
NBS  5 
NES  7 
NEND  7 


Miscellaneous  closed-surface  3-D  geometries  can  be  read  in  under  the  option  of  “3”  for 
the  variable  “IGEO”.  This  will  permit  easy  expansion  of  the  code  for  analyzing  other  closed- 
surface,  single-zone  PEC  geometries  other  than  a  sphere,  ogive,  and  cone-sphere.  Possible 
geometries  include  the  EMCC-defined  double  ogive,  cone-sphere  with  gap,  or  almond.  The  grid 
generated  for  the  geometry  must  incorporate  several  important  features.  First,  the  value  for  “KL” 
must  be  odd.  Several  sections  of  the  code  assume  that  KL  is  odd  to  correctly  generate  grid 
overlap  areas  and  to  fill  in  voids  in  the  calculation  of  the  fluxes  .  The  values  of  “IL”  and  “KL” 
can  be  odd  or  even.  Second,  the  code  assumes  that  the  PEC  surface  of  the  geometry  is  the  first 
cell  in  the  radial  direction  (spherical  coordinates).  The  PEC  surface  boundary  conditions  are 
imposed  at  the  first  cell  in  the  radial  direction.  Third,  the  FVTD  code  assumes  the  z  axis  is  the 
line  of  singularity,  or  rotation  (([)),  for  the  grid.  For  example,  the  tips  of  the  ogive  in  the  ogive 
grid  are  located  on  the  line  of  singularity.  The  geometry  doesn’t  necessarily  have  to  be  a  body- 
of-revolution  (BOR)  for  an  appropriate  grid  to  be  generated  correctly  for  the  code.  A  grid  could 
be  generated  for  the  NASA  almond  which  is  not  a  BOR  [55].  If  the  z  axis  is  not  the  line  of 
singularity,  the  grid  coordinates  must  be  rotated.  These  three  grid  requirements  are  necessary  for 
the  correct  use  of  the  FVTD  code. 

The  constitutive  parameters,  REP  and  RMU,  can  be  specified  only  for  real  values. 
Currently  the  code  cannot  incorporate  complex  values  for  the  permittivity  and  permeability 
although  this  would  be  a  relatively  simple  modification. 
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The  specification  of  the  Fourier  transform  variables  require  clarification.  For  use  with  a 
sinusoid  incident  wave  without  the  RCS  convergence  check,  the  “NBS”  and  “NFS”  variables 
specify  the  value  of  the  periods  in  which  the  Fourier  data  will  be  taken  which  is  usually  for  two 
to  three  periods  after  the  transients  diminish.  If  a  sinusoid  incident  wave  is  used  with  the  RCS 
convergence  check,  the  variables  “NBS”  and  “NES”  are  specified  in  the  same  way  except  the 
Fourier  transform  data  is  taken  for  only  one  period.  “NBS”can  begin  after  one  or  two  periods 
and  the  simulation  will  not  end  until  convergence  has  been  reached.  “NES”  and  “NEND”  can  be 
used  to  “override”  the  convergence  check  and  end  the  simulation.  For  a  Gaussian  pulse,  “NBS” 
is  not  used  since  the  Fourier  transform  data  is  taken  beginning  with  the  first  time  step.  “NES” 
and  “NEND”  are  specified  to  “override”  the  threshold  check.  The  simulation  ends  if  the 
amplitude  of  the  scattered  field  is  less  than  140  dB  below  the  peak  of  the  incident  field;  however, 
if  the  “NES”  value  is  reached  before  this  point,  the  simulation  ends.  The  override  gives  the  user 
a  method  to  end  the  test  and  obtain  data  without  running  out  of  computer  queue  time. 

If  a  Gaussian  pulse  is  used  for  the  incident  wave,  several  parameters  must  be  specified. 
Table  B.l  is  included  as  a  reference  for  specifying  the  Gaussian  pulse  parameters.  First,  the 
bandwidth  of  the  pulse  is  required.  Typically,  only  one-third  of  the  bandwidth  of  the  pulse  is 
used  [18].  Therefore,  a  pulse  must  be  specified  which  has  three  times  the  bandwidth  of  the 
desired  bandwidth.  For  example,  if  the  RCS  values  are  required  for  6  GHz,  the  parameters  for  a 
pulse  with  a  bandwidth  of  18  GHz  (f)  must  be  specified.  The  code  was  modified  to  incorporate 
up  to  15  different  frequencies;  therefore,  “WD”  must  be  specified  accordingly  as  to  not  exceed 
the  limit  of  15  frequencies  within  the  specified  bandwidth.  “GPD”  is  always  specified  as  4.0  so 
the  Gaussian  pulse  is  truncated  at  140  dB  down  from  the  peak.  “GPD”  refers  to  the  Gaussian 
pulse  delay  and  should  not  be  confused  with  GPD  used  in  Chapter  4.  The  transients  produced  at 
this  truncation  level  are  negligible.  The  value  of  to ,  computed  in  the  code,  will  be  the  value  as 
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shown  in  the  table..  The  value  of  “W”  should  be  specified  for  the  frequency  in  which  the  grid  is 
optimized.  The  values  in  the  table  are  normalized  to  the  velocity  of  a  wave  in  the  medium  and 


are  specified  in  the  same  way  in  the  input  file.  The  parameters  in  the  FVTD  code  are  also 
normalized  to  the  velocity  of  the  medium  (i.e.  velocity  of  an  electromagnetic  wave  in  free  space). 


Table  B.l:  Gaussian  Pulse  Parameters  for  a  Specified  Bandwidth 
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B.2.2  Ogive  Grid  Modifications 


The  majority  of  the  simulations  were  performed  for  the  ogive.  A  detailed  grid  resolution 
study  was  performed.  Addition  and  deletion  of  grid  points  increases  the  accuracy  and  the  ability 
to  obtain  grid  resolution  requirements  for  the  tips  of  the  ogive.  The  subroutine  for  the  ogive  grid 
is  listed  below.  A  similar  but  shorter  subroutine  was  used  for  the  cone-sphere  grid. 


SUBROUTINE  OGIVE (NG PD) 

PARAMETER { ID=75 , JD=125 , KD=45 , KF=15 , PI=3 . 14159265359 ) 

COMMON  /CO/  SXI(ID, JD,KD,3) ,SET(ID,JD,KD,3) ,SZT(ID, JD,KD,3) , 

1  V(ID, JD,KD) ,X(ID,JD,KD) ,Y(ID, JD,KD) ,Z(ID, JD,KD) 

COMMON  /IN/  AMPX,AMPY,AMPZ,DEL,GPD,GPT, ISRC,W,WG(KF) ,XDIS,yDIS, 
1  ZDIS,XK,YK, ZK 

COMMON  /PR/  AMP,CFL,DT,EPS,IB,IL,JL,KL.ILM, JLM,KLM,KWM,REP,RMU, 
1  RK,RO,RP,SS,T 

DIMENSION  XG(ID, JD,KD) , YG { ID, JD, KD) , ZG ( ID, JD, KD) 

CHARACTER  *72  HEADER 
C 

C  VARIABLES  USED  ONLY  IN  THE  OGIVE  SUBROUTINE: 

C  DJ:  DISTANCE  IN  THE  J  DIRECTION 

C  DPH:  ANGLE  INCREMENT  IN  THE  PHI  DIRECTION 

C  GPD'S:  GRID  POINT  DENSITIES  IN  EACH  DIRECTION 

C  I,J,KG:  SIZE  OF  ORIGINAL  GRID 

C  IRI:  FACTOR  TO  INCREASE  NUMBER  OF  POINTS  IN  RADIAL  DIRECTION 
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o  o  o  o  o  on 


C  I  AND  J  INDICES  DURING  AN  INTERMEDIATE  CALCULATION 

C  JTD:  THETA  INDEX  TO  WHICH  TIP  DENSITY  WILL  BE  INCREASED 

C  JTHD:  FACTOR  TO  DECREASE  NUMBER  OF  POINTS  IN  THETA  DIRECTION 

C  JTHI:  FACTOR  TO  INCREASE  NUMBER  OF  POINTS  IN  THETA  DIRECTION 

C  NGPD:  NUMBER  OF  POINTS  REQUIRED  FOR  OGIVE  TIPS 

C  NR:  EXTRA  NUMBER  OF  GRID  POINTS  ADDED  IN  RADIAL  DIRECTION 

C  R,D:  DISTANCE  TO  A  GRID  POINT  OR  BETWEEN  GRID  POINTS 

C  ROE:  EXTENSION  OF  OUTER  BOUNDARY  IN  THE  RADIAL  DIRECTION 

C  WL :  WAVELENGTH 

C  X,Y,ZG:  GRID  POINTS  AT  AN  INTERMEDIATE  CALCULATION 

C 

C  SPECIFY  VARIOUS  FORMATS  FOR  FILES 

100  FORMAT (E13 . 6 , IX, E13 . 6 , IX, E13 . 6 ) 

101  FORMAT (A72) 

102  FORMAT (Ell. 4, IX, Ell. 4, IX, Ell. 4) 

C 

IG=10 
JG=121 
KG=3  0 

WL=2.0*PI/W 
RO=RO*WL 
ROE=RO-0 . 08 
JL4=JL-4 
KL4=KL-4 

DPH=2 .0*PI/FLOAT(KL-5) 

INPUT  THE  OGIVE  GRID  (UNITS  ARE  IN  METERS) 

OPEN (UNIT=2 , FILE= ' ogive . grd ' , STATUS= ' OLD ' ) 

READ (2, 101)  HEADER 
READ (2, 101)  HEADER 
READ (2, 101)  HEADER 
DO  1  K=1,KG 
DO  1  J=1,JG 
DO  1  1=1, IG 

READ (2, 100)  X(I, J,K) ,Y(I, J,K) ,Z(I, J,K) 

1  CONTINUE 
CLOSE (2) 

ADJUST  GRID  TO  DESIRED  DIMENSIONS 
CHANGE  NUMBER  OF  PHI  GRID  POINTS  (K) 

DO  2  1=1, IG 
DO  2  J=1,JG 
R=0.0 
D=0 . 0 

DO  2  K=1,KG 

R=R+SQRT(X(I, J,K) **2+Y(I, J,K) **2) 

D=D+Z (I, J,K) 

IF(K.EQ.KG)  THEN 
R=R/ FLOAT (KG) 

D=D/ FLOAT (KG) 

DO  3  KOG=l,KL4 

PH=FLOAT (KOG-1) *DPH 
XG ( I , J, KOG) =R*COS (PH) 

YG(I, J,KOG) =R*SIN(PH) 

ZG(I, J,KOG) =D 
CONTINUE 
END  IF 
CONTINUE 

CHANGE  NUMBER  OF  R  GRID  POINTS  (I) 
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NR=NINT{ROE/ { ZG ( IG, 1 , 1 ) -ZG { (IG-1) ,1,1) ) ) 

IF{NR.GT.O)  THEN 
DO  5  K=1,KL4 
DO  5  J=1,JG 

DX=XG{IG, J,K) -XG( (IG-l) , J,K) 

DY=YG(IG, J,K) -YG( (IG-l) , J,K) 

DZ=ZG(IG, J,K) -ZG( (IG-l) , J,K) 

DO  5  1=1, NR 

XG( (IG+I) , J,K)=FLOAT(I) *DX+XG ( IG, J, K) 

YG( (IG+I) , J,K) =FLOAT(I) *DY+YG ( IG, J, K) 

ZG( (IG+I) , J,K) =FLOAT(I) *DZ+ZG ( IG, J, K) 

5  CONTINUE 
END  IF 

C 

IRI=NINT (FLOAT (ID) /FLOAT ( IG+NR) ) 

IF(IRI.EQ.O)  IRI=1 
IL= ( (IG-l) +NR) *IRI+1 
DO  6  K=1,KL4 
DO  6  J=1,JG 
DO  6  1=1, IL 

Cl=FLOAT(MOD( (I-l) ,IRI) ) /FLOAT (IRI) 

11= (I-l) /IRI 

X(I, J,K) =C1* (XG( (11+2) , J,K) -XG( (Il+l) , J,K) ) +XG ( (Il+l) , J,K) 
Y(I, J,K) =C1* (YG( (11+2) , J,K) -YG( (Il+l) , J,K) ) +YG ( (Il+l) , J,K) 
Z(I, J,K) =C1* (ZG( (11+2) , J,K) -ZG( (Il+l)  ,  J,K) ) +ZG( (Il+l) , J,K) 

6  CONTINUE 
C 

C  CHANGE  NUMBER  OF  THETA  GRID  POINTS  (J) 

IF(JL4.LT. JG)  THEN 

JTHD=NINT ( FLOAT ( JG ) / FLOAT ( JL4 ) ) 

JTHI=0 

JL=( (JG-18) /JTHD)+23 
JL4=JL-4 
DO  10  K=1,KL4 
DO  10  J=1,JL4 
DO  10  1=1, IL 

IF(J.LE.IO)  THEN 
XG(I, J,K)=X(I, J,K) 

YG(I, J,K)=Y{I, J,K) 

ZG(I, J,K)=Z(I,J,K) 

ELSE  IF(J.GE. (JL4-9) )  THEN 
J1=JG- (JL4-J) 

XG(I, J,K) =X(I, J1,K) 

YG(I, J,K) =Y(I, J1,K) 

ZG(I, J,K)=Z(I, J1,K) 

ELSE 

XG(I, J,K) =X(I, ( (J-9) *JTHD+9- (JTHD-1) ) , K) 

YG(I, J,K) =Y(I, ( (J-9) *JTHD+9- (JTHD-1) ) ,K) 

ZG(I, J,K) =Z (I, { (J-9) *JTHD+9- (JTHD-1) ) ,K) 

END  IF 

10  CONTINUE 
END  IF 
C 

IF(JL4.EQ. JG)  THEN 
DO  11  K=1,KL4 
DO  11  J=1,JL4 
DO  11  1=1, IL 

XG(I, J,K)=X(I, J,K) 

YG(I, J,K)=Y(I, J,K) 

ZG(I, J,K)=Z(I,J,K) 
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11  CONTINUE 
END  IF 

IF{JL4.GT. JG)  THEN 
JTHD=0 

JTHI=NINT {FLOAT (JL4} /FLOAT (JG) ) 

JL=(JG*JTHI) - (JTHI-l)+4 
JL4=JL-4 
DO  12  K=1,KL4 
DO  12  J=1,JL4 

Cl=FLOAT (MOD ( (J-1) , JTHI) ) /FLOAT(JTHI) 

J1=(J-1) /JTHI 
DO  12  1=1, IL 

XG(I, J,K) =C1* (X(I, {Jl+2) ,K) -X{I, (Jl+1) ,K) ) +X(I, (Jl+1) , K) 
YG(I, J,K) =C1* (Y(I, (Jl+2) ,K) -Y(I, (Jl+1) ,K) ) +Y(I, (Jl+1) ,K) 
ZG{I, J,K) =C1* (Z (I, (Jl+2) ,K)-Z(I, (Jl+1) ,K) )+Z{I, (Jl+l) ,K) 

12  CONTINUE 
END  IF 

INCREASE  GRID  POINT  DENSITY  AT  TIPS  OF  OGIVE 
JTD=10 

D=SQRT{ (XG{1,2, 1) -XG{1,1,1) ) **2+ ( YG ( 1 , 2 , 1 ) -YG { 1 , 1 , 1 ) ) **2 
1  +{ZG{1,2,1)-ZG(1,1,1) )**2) 

NGPD=NINT { FLOAT (NGPD) *D/WL) 

IF(NGPD.EQ. 0)  NGPD=1 
JL=JL+2* (NGPD-1) * (JTD-l) 

JL4=JL-4 

DO  16  J=l, { (JTD-l) *NGPD+1) 

Cl=FLOAT{MOD{ (J-1) ,NGPD) ) /FLOAT (NGPD) 

J1=(J-1) /NGPD 
DO  16  K=1,KL 
DO  16  1=1, IL 

X(I, J,K) =C1* (XG(I, (Jl+2) ,K)-XG(I, (Jl+1) ,K) ) +XG ( I , (Jl+1) ,K) 
Y(I, J,K) =C1* {YG(I, (Jl+2) ,K) -YG(I, (Jl+1) ,K) ) +YG ( I , (Jl+1) ,K) 
Z (I, J,K) =C1* {ZG(I, {Jl+2) ,K) -ZG(I, (Jl+1) ,K) ) +ZG(I, (Jl+l) ,K) 
;  CONTINUE 

DO  17  J= ( (JTD-l) *NGPD+2) , (JL4- ( (JTD-l) *NGPD+1) ) 

J1=J- (JTD-l) * (NGPD-l) 

DO  17  K=1,KL 
DO  17  1=1, IL 

X(I, J,K)=XG(I, J1,K) 

Y(I, J,K) =YG(I, J1,K) 

Z{I, J,K)=ZG(I,J1,K) 

CONTINUE 

J2=JL4- ( ( JTD-l) *NGPD) 

DO  18  J=J2,JL4 

Cl=FLOAT{MOD( (J-J2) ,NGPD) ) /FLOAT (NGPD) 

Jl= ( J-J2) /NGPD+ ( J2- (JTD-l) * (NGPD-1) ) 

DO  18  K=1,KL 
DO  18  1=1, IL 

X(I, J,K)=C1* {XG(I, (Jl+l) ,K) -XG{I, J1,K) )+XG{I,Jl,K) 

Y(I, J,K) =C1* {YG(I, (Jl+1) ,K) -YG(I, J1,K) ) +YG(I, J1,K) 

Z (I, J,K) =C1* (ZG(I, (Jl+1) ,K) -ZG(I, J1,K) )+ZG(I,Jl,K) 

CONTINUE 

OUTPUT  THE  GRID  FOR  PLOTTING 

OPEN {UNIT=3,FILE=' ogive. pit' , STATUS =' UNKNOWN ' ) 
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WRITE (3,*)  ' TITLE=OGIVE ' 

WRITE (3,*)  ' VARIABLES=X,Y,Z' 

WRITE(3,*)  'ZONE  T=" GRID",  I=',IL,',  J=',JL4,',  K=',KL4 
DO  20  K=1,KL4 
DO  20  J=1,JL4 
DO  20  1=1, IL 

WRITE (3, 102)  X(I, J,K) ,Y(I, J,K) ,Z(I, J,K) 

20  CONTINUE 
CLOSE (3) 


CALCULATE  GRID  POINT  DENSITIES 
J=JL/2+l 
KL5=KL-5 

GPDKB=KL5/ (W*X(1, J, 1) ) 

GPDKT=KL5 / ( W*X (1,2,1)) 

GPDII=WL/ (X{2, (JL/2+1) ,1)-X(1, (JL/2+1) ,1) ) 
GPDIO=WL/ (X(IL, (JL/2+1) ,1)-X( (IL-1) , (JL/2+1) ,1) ) 


D=0 . 0 

DO  30  J= ( ( JTD-1) *NGPD+2) , (JL4- (JTD-1) *NGPD) 

D=SQRT( (X(l, J, 1) -X(l, (J-1) ,1) ) **2+(Z(l, J, 1) -Z (1, (J-1) , 1) ) **2) +D 
3  0  CONTINUE 

GPD JB=WL* FLOAT ( (JL4- (JTD-1 ) *NGPD) - ( (JTD-1) *NGPD+2 ) ) /D 


D=0.0 

DO  31  J=2, ( (JTD-1) *NGPD+1) 

D=SQRT( (X(l, J, 1) -X(l, (J-1) ,1) ) **2+(Z(l, J, 1) -Z (1, (J-1) , 1) ) **2) +D 
CONTINUE 

GPD JT=WL * FLOAT ( (JTD-1) *NGPD-1) /AMAX1(D,SS) 


CALCULATE  SURFACE  LENGTH  IN  J  DIRECTION 
DJ=0.0 

DO  32  J=2,JL4 
J1=J-1 

DJ=SQRT( (X(l, J, 1) -X(l, Jl,l) ) **2+(Z(l, J, 1) -Z (1, Jl, 1) ) **2) +DJ 
32  CONTINUE 
C 

PRINT  *, 'OGIVE:  RO , ROE , NR, IRI , IL  ' , RO, ROE , NR, IRI , IL 

PRINT  * , ' OGIVE :  JTHD , JTHI , JTD , NGPD , JL  ' , JTHD , JTHI , JTD , NGPD , JL 

PRINT  *, 'OGIVE:  GPD=KB , KT, II , lO, JB, JT  ’,  GPDKB , GPDKT , GPDII , 

1  GPDIO,GPDJB,GPDJT 

PRINT  *,'L,RX,DJ  ' ,Z(1,1,1)-Z(1,JL4,1) ,X(1, (JL/2+1) ,1) ,DJ 
C 

RETURN 

END 

Q  ********************************************************************* 


B.2.3  Incident  Field 

The  initial  parameters  for  the  incident  field  are  calculated  prior  to  the  time  loop.  The 
components  of  the  E  field  vector,  displacement  vector,  and  the  propagation  unit  vector  are 
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calculated  from  the  input  parameters.  The  frequencies  are  also  calculated  if  the  Gaussian  pulse 
incident  wave  is  specified.  If  a  Gaussian  pulse  is  used,  a  delay  is  introduced  so  the  truncated 
edge  of  the  pulse  is  initially  located  at  the  point  on  the  surface  of  the  object  farthest  away  from 
the  origin  in  the  direction  of  the  incident  wave.  The  origin  is  located  inside  of  the  PEC  scatterer. 
The  section  of  the  code  which  calculates  the  initial  incident  field  parameters  is  shown  below. 


C  CALCULATE  INCIDENT  FIELD  PARAMETERS 
C  CONVERT  ANGLES  FROM  DEGREES  TO  RADIANS 
COSPH=COS (PHI*PI/180.0) 

SINPH=SIN(PHI*PI/180 .0) 

COSTH=COS ( THI *PI/180.0) 

SINTH=SIN(THI*PI/180 . 0) 

FIND  THE  COMPONENTS  OF  THE  E  FIELD  VECTOR 
AMPX=AMP* (ETHI*COSTH*COSPH-EPHI*SINPH) 

AMPY=AMP* (ETHI*COSTH*SINPH+EPHI*COSPH) 

AMPZ=AMP* (-ETHI*SINTH) 

FIND  THE  COMPONENTS  OF  THE  DISPLACEMENT  UNIT  VECTOR  FOR  DELAY 
XDIS=SINTH*COSPH 
YDIS=SINTH*SINPH 
ZDIS=COSTH 

CALCULATE  THE  ANGLES  AND  COMPONENTS  OF  K  (PROPAGATION)  UNIT  VECTOR 
PHK=PHI-180 . 0 
THK=-THI+180 . 0 

XK=SIN(THK*PI/180 . 0) *COS (PHK*PI/180 . 0) 

YK=SIN(THK*PI/180 .0) *SIN (PHK*PI/180 . 0 ) 

ZK=COS ( THK*  PI / 180.0) 

INITIALIZE  THE  FREQUENCIES 
DO  4  KW=1,KF 
WG(KW)=0.0 
CONTINUE 

ASSIGN  THE  FREQUENCY  OF  THE  SINUSOID 
KVM=1 
WG ( 1 ) =W 

CALCULATE  THE  PARAMETERS /FREQUENCIES  FOR  THE  GAUSSIAN  PULSE 
IF  (ISRC.EQ.l)  THEN 
GTH=1 . OE-8 

WM=2*PI* (1 . 0/3 . 0) * (1/ (PI*GPT) ) *SQRT(-LOG(1.0E-7) ) 

KWM=INT(1 . 0/WD*WM) +1 
PRINT  *, 'W,WD,KWM  ARE:  ',W,WD,KWM 
DO  5  KW=1,KWM 

WG (KW) =FLOAT (KW-1) *WD 
CONTINUE 
END  IF 

CALCULATE  DELAY  SO  INCIDENT  WAVE,  IF  GAUSSIAN,  PROPAGATES  CORRECTLY 
DEL=0 . 0 
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IF  (ISRC.EQ.l)  THEN 
DO  6  K=1,KLM 
DO  6  J=1,JLM 

D=XC1 ( J,K) *XDlS+yCl (J,K) *YDIS+ZC1 (J,K) *ZDIS 
IF  (D.GT.0.0)  THEN 
DEL=AMAX1 (DEL,D) 

END  IF 

6  CONTINUE 

END  IF 


The  values  of  the  incident  field,  either  from  a  sinusoid  or  Gaussian  pulse  incident  wave, 
at  the  surface  of  the  scatterer  are  calculated  in  the  subroutine  SOURCE.  The  code  listing  for  the 
SOURCE  subroutine  is  shown  below. 


c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


1 


1 

1 

1 


SUBROUTINE  SOURCE (BIXl , BIX2 , BIYl , BIY2 , BIZI , BIZ2 , DIXl , DIX2 , DIYl , 
DIY2,DIZ1,DIZ2, J,K,TF) 

PARAMETER ( ID=75 , JD=125 , KD=45 , KF=15 ) 

COMMON  /CC/  XCl ( JD,KD) ,XC2 (JD,KD) , YCl (JD,KD) , YC2 ( JD,KD) , 
ZC1(JD,KD) ,ZC2(JD,KD) 

COMMON  /IN/  AMPX,AMPY,AMPZ,DEL,GPD,GPT, ISRC,W,WG(KF) ,XDIS,YDIS, 
ZDIS,XK,YK, ZK 

COMMON  /PR/  AMP,CFL,DT,EPS,IB,IL, JL,KL,ILM, JLM,KLM,KWM,REP,RMU, 
RK,RO,RP,SS,T 


VARIABLES 
BIXl, 2 : 
BIYl, 2: 
BIZI, 2 ; 
DIS1,2 : 
DIXl, 2 : 
DIYl, 2 : 
DIZl, 2 : 
GAUS1,2: 
SINWl, 2 : 


USED  ONLY  IN  THE  SOURCE  SUBROUTINE: 

X  COMPONENTS  OF  INCIDENT  B  FIELD  AT  CELL  CENTERS  1,2 

Y  COMPONENTS  OF  INCIDENT  B  FIELD  AT  CELL  CENTERS  1,2 

Z  COMPONENTS  OF  INCIDENT  B  FIELD  AT  CELL  CENTERS  1,2 

SPATIAL  DELAY  AT  A  CELL  WITH  CENTER  OF  1=1,2 

X  COMPONENTS  OF  INCIDENT  B  FIELD  AT  CELL  CENTERS  1,2 

Y  COMPONENTS  OF  INCIDENT  B  FIELD  AT  CELL  CENTERS  1,2 

Z  COMPONENTS  OF  INCIDENT  B  FIELD  AT  CELL  CENTERS  1,2 

INCIDENT  FIELD  AS  A  FCT  OF  W, GPT , GPD, DEL ,  AND  DISl , 2 
INCIDENT  FIELD  AS  A  FCT  OF  W,T,  AND  DISl, 2 


C1=SQRT(RMU/REP) 


CALCULATE  SPATIAL  DELAY 

DIS1=XC1 ( J, K) *XDIS+YC1 ( J, K) *YDIS+ZC1 ( J, K) *ZDIS 
DIS2=XC2 ( J, K) *XDIS+YC2 { J, K) *YDIS+ZC2 ( J, K) *ZDIS 
C 

C  CALCULATE  D  FIELD  IF  A  SINUSOID  IS  REQUIRED 
IF  (ISRC.EQ.O)  THEN 

SINW1=SIN{W* (TF+DISl) ) 

SINW2=SIN(W* (TF+DIS2) ) 

DIX1=REP*AMPX*SINW1 
DIX2=REP*AMPX*SINW2 
DIY1=REP*AMPY* SINWl 
DIY2=REP*AMPY*SINW2 
DIZ1=REP*AMPZ*SINW1 
DIZ2=REP*AMPZ*SINW2 

C 

C  CALCULATE  D  FIELD  IF  A  GAUSSIAN  PULSE  IS  REQUIRED 
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ELSE 

RGPT=1.0/GPT 

GP=GPD*GPT 

GAUS1=EXP (- ( (TF+DISl-DEL-GP) *RGPT) **2) 

GAUS2=EXP(-( (TF+DIS2-DEL-GP) *RGPT) **2) 

DIX1=REP*AMPX*GAUS1 
DIX2=REP*AMPX*GAUS2 
DIY1=REP*AMPY*GAUS1 
DIY2=REP*AMPY*GAUS2 
DIZ1=REP*AMPZ*GAUS1 
DIZ2=REP*AMPZ*GAUS2 
END  IF 
C 

C  CALCULATE  THE  COMPONENTS  OF  THE  B  FIELD 
BIX1=C1* {YK*DIZ1-ZK*DIY1) 

BIX2=C1* (YK*DIZ2-ZK*DIY2) 

BIY1=C1* (ZK*DIX1-XK*DIZ1) 

BIY2=C1* (ZK*DIX2-XK*DIZ2) 

BIZ1=C1* (XK*DIY1-YK*DIX1) 

BIZ2=C1* {XK*DIY2-YK*DIX2) 

C 

RETURN 

END 

********************************************************************* 


As  can  be  seen  in  the  code,  the  spatial  and  time  delay  is  taken  into  account  for  the 
sinusoid  and  Gaussian  pulse  incident  wave.  The  location  of  the  grid  point  on  the  surface  of  the 
scatterer  is  accounted  for  in  the  calculation  of  the  sinusoid  incident  wave.  For  the  Gaussian 
pulse,  the  time  delay  for  the  center  of  the  pulse,  the  time  delay  due  to  the  length  of  the  object, 
and  the  time  delay  due  to  the  location  of  the  grid  point  on  the  surface  of  the  scatterer  are  used  in 
calculating  the  amplitude  of  the  pulse. 


B.2.4  RCS  Convergence  Check 

The  RCS  convergence  check  is  used  to  end  the  time  loop  when  the  RCS  values  have 
converged  to  within  a  certain  value  when  using  a  sinusoid  incident  wave.  The  RCS  values  are 
sampled  every  10  degrees  every  period.  If  the  values  are  with  0.1  dB  of  the  previously  sampled 
values,  convergence  has  been  reached  and  the  time  loop  ends.  The  frequency  data  for  the 
scattered  fields  are  set  to  zero  prior  to  each  period. 
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Q  **************************************************************************** 

SUBROUTINE  RCSCONV ( FS , FU , N, NEND, NES , NP , RCSEC , RI ) 
PARAMETER(ID=75,JD=125,KD=45,PI=3 . 141592 653 59, KF=15) 

COMMON  /CC/  XCl (JD,KD} ,XC2 (JD,KD) ,YC1(JD,KD) ,YC2 (JD,KD) , 

1  ZCl (JD,KD) , ZC2 (JD,KD) 

COMMON  /CO/  SXKID,  JD,KD,3)  ,SET(ID,  JD,KD,3)  ,SZT(ID,  JD,KD,3)  , 

1  V(ID, JD,KD) ,X(ID, JD,KD) ,Y(ID, JD,KD) , Z (ID, JD,KD) 

COMMON  /IN/  AMPX,AMPY,AMPZ,DEL,GPD,GPT, ISRC,W,WG(KF) ,XDIS,YDIS, 

1  ZDIS,XK, YK, ZK 

COMMON  /PR/  AMP,CFL,DT,EPS,IB,IL, JL,KL,ILM, JLM,KLM,KWM,REP,RMU, 

1  RK,RO,RP, SS,T 

COMPLEX  AX,AY,AZ,CX,CY,CZ,DE,FS(JD,KD,6,KF) , FU ( JD, KD, 6 , KF) , 

1  FUEC ( JD , KD , 6 ) , RI , SUMEKl , SUMEK2 , SUMEK3 , SUME JKl , SUME JK2 , 

2  SUMEJK3 
DIMENSION  RCSEC (19) 

C 

C  VARIABLES  USED  ONLY  IN  THE  RCSCONV  SUBROUTINE: 

C  SEE  MAIN  PROGRAM  AND  RCS  SUBROUTINE  FOR  MOST  VARIABLE  DESCRIPTIONS 

C  FUEC:  FREQUENCY  DATA  USED  FOR  CONVERGENCE  CHECK 

C  RCSCH:  RATIO  OF  CURRENT  RCS  VALUES  TO  PREVIOUS  RCS  VALUES 

C  RCSCHI:  INTERMEDIATE  VALUE  FOR  RCSCH 
C 

SME=SQRT (RMU/REP) 

Cl=l. 0/REP 
C2=1.0/RMU 

CALCULATE  INTEGRATION  LIMITS 
JLM1=JLM-1 
KLM4=KLM-4 

NORMALIZE  SCATTERED  FIELD  DATA  TO  INCIDENT  FIELD 
DO  5  K=1,KL 
DO  5  J=1,JL 

RFSB=SQRT(CABS (FS (J,K, 1, 1) ) **2+CABS(FS (J,K,2,1))**2 
1  +CABS(FS(J,K,3,1) )**2) 

DO  6  L=l,3 

FUEC ( J,K,L) =FU(J,K,L, 1) /RFSB 

6  CONTINUE 

RFSD=SQRT(CABS (FS (J,K, 4, 1) ) **2+CABS (FS (J,K,5,1))**2 
1  +CABS(FS(J,K,6,1))**2) 

DO  7  L=4,6 

FUEC(J,K,L)=FU(J,K,L,1) /RFSD 

7  CONTINUE 
5  CONTINUE 

GENERATE  THE  BISTATIC  RADAR  CROSS  SECTION 


PHO=0 . 0 
SINPO=SIN(PHO) 

COSPO=COS (PHO) 

CHECK  FOR  CONVERGENCE  OF  RCS  VALUES 
RCSCH=0 . 0 
DO  10  JO=l,19 

THO= (PI*FLOAT( JO-1) *10.0) /180.0 
SINTO=SIN(THO) 

COSTO=COS (THO) 
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CALCULATE  CARTESIAN  COMPONENTS  OF  THE  R  (OBSERVATION)  UNIT  VECTOR 
BX=SINTO*COSPO 
BY=SINTO*SINPO 
BZ=COSTO 

PERFORM  THE  SURFACE  INTEGAL  (FREE  SPACE  GREEN'S  FUNCTION) 

SET  INITIAL  VALUES  FOR  THE  SURFACE  INTEGRAL, J  &  K  LOOPS 
SUMEJK1= (0 . 0, 0 .0) 

SUMEJK2=(0.0,0.0) 

SUMEJK3=(0.0,0.0) 

DO  20  J=3,JLM1 

SET  INITIAL  VALUES  FOR  THE  SURFACE  INTEGRAL,  K  LOOP 
SUMEK1= (0 . 0, 0 .0) 

SUMEK2= (0 . 0, 0 .0) 

SUMEK3=(0.0,0.0) 

DO  30  K=1,KLM4 

CALCULATE  THE  DOT  PRODUCT  OF  R'  AND  R 

PV=BX*XC2 ( J,K) +BY*YC2 (J,K) +BZ*ZC2 (J,K) 

WP=WG(1) *PV 

GENERATE  SURFACE  INTEGRAL  SCALE  &  METRICS  OF  OUTWARD  NORMAL 
DS=SQRT( (0.5* (SXI (IB, J,K, 1) +SXI ( IB-1 , J, K, 1 ) ) ) **2 

1  +(0.5* (SXI(IB, J,K,2)+SXI(IB-1, J,K,2) ) ) **2 

2  +(0.5* (SXI(IB, J,K,3)+SXI(IB-1, J,K,3) ) ) **2) 
SUMX=AMAX1(DS,SS) 

RSUMX=1.0/SUMX 
XIX=SXI (IB, J,K, 1) *RSUMX 
XIY=SXI (IB, J,K,2) *RSUMX 
XIZ=SXI(IB, J,K,3)*RSUMX 

GENERATE  THE  CROSS  PRODUCTS  OF  THE  ELECTROMAGNETIC  FIELD 
THE  SURFACE  ELECTRIC  &  MAGNETIC  CURRENT  DENSITIES 
FOR  THE  EQUIVALENCE  THEOREM 

AX=C1* (XIY*FUEC(J,K, 6) -XIZ*FUEC(J,K, 5) ) 

AY=C1* (XIZ*FUEC (J,K, 4) -XIX*FUEC (J,K, 6) ) 

AZ=C1* (XIX*FUEC(J,K, 5) -XIY*FUEC(J,K,  4)  ) 

CX=C2* (XIY*FUEC(J,K,3)-XIZ*FUEC(J,K,2) ) 

CY=C2* (XIZ*FUEC(J,K,1)-XIX*FUEC(J,K,3) ) 

CZ=C2* (XIX*FUEC(J,K,2) -XIY*FUEC(J,K, 1) ) 

DE=C1* (XIX*FUEC ( J, K, 4) +XIY*FUEC ( J, K, 5 ) +XIZ*FUEC ( J, K, 6 ) ) 

IF( (J.EQ.3) .OR. (J.EQ.JLMl) )  THEN 
AX=AX/2.0 
AY=AY/2 . 0 
AZ=AZ/2 . 0 
CX=CX/2 . 0 
CY=CY/2 . 0 
CZ=CZ/2 . 0 
DE=DE/2 . 0 
END  IF 

PERFORM  THE  SURFACE  INTEGRATIONS  FOR  NEAR  TO  FAR  FIELD  TRANSFORM 
SUMEK1=SUMEK1+WG(1) * (DE*BX+ (AY*BZ-AZ*BY) -SME*CX) 

1  *CEXP(- (RI*WP) ) *DS 

SUMEK2=SUMEK2+WG (1) * (DE*BY+ (AZ*BX-AX*BZ ) -SME*CY) 
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*CEXP{- (RI*WP) ) *DS 

SUMEK3=SUMEK3+WG (1) * (DE*BZ+ (AX*BY-AY*BX) -SME*CZ) 

1  *CEXP{- (RI*WP) ) *DS 

C 

3  0  CONTINUE 

C 

SUME JKl  =  SUME JKl  +  SUMEKl 
SUME JK2 = SUME JK2  +  SUMEK2 
SUME JK3  =  SUME JK3  +  SUMEK3 
C 

2  0  CONTINUE 

C 

C  CALCULATE  THE  RCS:  4*PI * (R* *2 ) * | Es | * *2 / | Ei | * *2 
C  NOTE:  NORMALIZATION  WRT  THE  INCIDENT  FIELD  WAS  DONE  DURING  FT 
RCSE=RCSEC(JO) 

RCSEC (JO)  =  (CABS (SUME JKl) * *2+CABS ( SUMEJK2 ) *  * 2 +CABS ( SUME JK3 ) **2) 
1  /(4.0*PI) 

CHECK  FOR  CONVERGENCE,  CONTINUE  OR  END  TIME  LOOP 
RCSE=AMAX1 (ROSE, SS) 

RCSCHI=RCSEC (JO) /ROSE 
IF (RCSCHI . LT . 1 . 0 )  RCSCHI=1 . 0/RCSCHI 
RCSCH=AMAX1 (RCSCHI , RCSCH) 

10  CONTINUE 
C 

C  IF  ALL  RCS  VALUES  ARE  WITHIN  0.1  DB  OF  PREVIOUS  VALUES, 

C  CONVERGENCE  HAS  BEEN  REACHED 
IF (RCSCH.lt. 1.0233)  THEN 
NEND=N 
NES=NEND 

C  RESET  THE  INITIAL  VALUES  FOR  ALL  FOURIER  TRANSFORM  VARIABLES 
ELSE 

DO  40  L=l,6 
DO  40  K=1,KL 
DO  40  J=1,JL 

FU(J,K,L,1)=(0. 0,0.0) 

FS(J,K,L,1)=(0. 0,0.0) 

40  CONTINUE 

END  IF 

PRINT  *,'N, RCSCH:  ',N, RCSCH 
C 

RETURN 

END 


B.2.5  Scattered-Field  Threshold  Check 


The  time  loop  ends  if  the  amplitude  of  the  scattered  field  is  below  a  predetermined 
threshold  when  a  Gaussian  pulse  is  used  for  the  incident  wave.  The  amplitude  of  the  scattered 
field  on  the  virtual  surface  (I=IB)  is  checked  every  time  step  to  determine  the  maximum 
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amplitude  of  the  scattered  field.  The  time  loop  which  contains  the  threshold  check  is  listed 


below. 


C  PERFORM  THE  SOLVING  SEQUENCE 
N=1 

DO  WHILE (N.LE.NEND) 

T=FLOAT(N-l) *DT 

MR=0 

M=1 

^  •k'k-k-k'k'k-k-k-k'k'k'k-k 

CALL  BC(M,MR) 

DO  11  MR=1,4 
************ 

CALL  FXI 
CALL  GETA 
CALL  HZETA 
CALL  SUM (MR) 

★  ★★★★★★★★★★★ 

11  CONTINUE 

IF( (N.GT.NES) .OR. (N.LT.NBS) )  GOTO  10 
PERFORM  FOURIER  TRANSFORM 
TF=T+DT 
RU=0.0 

DO  12  KW=1,KWM 
WT=WG(KW) * (TF) 

DO  12  K=1,KL 
DO  12  J=1,JL 

********************************************************* 
CALL  SOURCE (BIX1,SN(1) ,BIY1,SN(2) ,BIZ1,SN(3) ,DIX1,SN(4) , 

1  DIY1,SN(5) ,DIZ1,SN(6) , J,K,TF) 

********************************************************* 

IF  GAUSSIAN  SOURCE,  CHECK  AMPLITUDE  OF  SCATTERED  FIELD,  END  TIME  LOOP 
IF  BELOW  THRESHOLD 

IF(ISRC.EQ.l)  THEN 

RUL=SQRT(U(1,IB, J,K,4) * *2+U ( 1 , IB , J , K, 5 ) **2 
1  +U(1,IB,J,K,6)**2) 

RU=AMAX1 {RU,RUL) 

END  IF 

DO  12  L=l,6 

FU(J,K,L,KW)=FU(J,K,L,KW) +(U(1, IB, J,K,L) +SN(L) ) 

1  *CEXP{RI*WT) 

FS (J,K,L,KW)=FS(J,K,L,KW)+SN(L) *CEXP{RI*WT) 

L2  CONTINUE 

IF{ (ISRC.EQ.l) .AND. (N.GE.lOO) )  THEN 
IF{RU.LE.GTH)  NEND=N 
NES=NEND 
END  IF 

IF  SINUSOID  SOURCE,  RUN  CONVERGENCE  CHECK  ON  RCS  VALUES 
IF ( ( ISRC . EQ . 0 ) . AND . ( ICON . EQ . 1 ) )  THEN 
IF {MOD(N,NP) .EQ. 0)  THEN 

Q  ********************************************** 
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c 


CALL  RCSCONV ( FS , FU , N, NEND , NES , NP , RCSEC , RI ) 

END  IF 
END  IF 
C 

10  N=N+1 

END  DO 

C  END  OF  TIME  LOOP 


B.2.6  Bistatic-to-Monostatic  Approximation 


A  description  of  the  bistatic-to-monostatic  approximation  is  given  in  Chapter  3.  The 
code  is  not  a  subroutine  for  the  FVTD  code  but  is  an  independent  program.  The  RCS  values 
from  the  “rcsmO.plt”  files  are  used  for  the  HH  polarization  approximation.  The  incident  wave 
must  be  specified  with  ETH  polarization  to  obtain  HH  polarization  results.  To  obtain  RCS 
results  for  the  VV  polarization,  the  “rcseO.plt”  files  are  used.  The  incident  wave  must  be 
specified  with  EPH  polarization  to  obtain  VV  polarization  results.  The  code  listed  below  is  for 
obtaining  a  monostatic  sweep  for  the  ogive  for  HH  polarization.  A  similar  file  is  used  for  VV 
polarization  and  to  obtain  the  monostatic  results  for  the  cone-sphere. 


C  THIS  FORTRAN  CODE  USES  CONTINUOUS  BISTATIC  VALUES  AND  SEVERAL 
C  MONOSTATIC  VALUES  TO  APPROXIMATE  A  CONTINUOUS  MONOSTATIC 
C  SOLUTION  FOR  THE  OGIVE 
C 

C  SOURCE:  THE  MONOSTATIC/BISTATIC  APPROXIMATION 
C  MICHAEL  SCHUH,  ALEX  WOO,  MICHAEL  SIMON 

C  IEEE  ANTENNAS  AND  PROPAGATION  MAGAZINE 

C  VOL.  36,  NO.  4,  AUGUST  1994,  PAGES  76-78 

C 

C  VARIABLE  DESCRIPTIONS 

C  DEG:  ANGLE  AT  WHICH  MONOSTATIC  VALUE  IS  APPROXIMATED 

C  NB:  INDEX  FOR  BISTATIC  RCS  VALUE 

C  NU:  INDEX  FOR  FILE  UNIT  WHICH  IS  CURRENTLY  BEING  USED 

C  RCSE:  MONOSTATIC  AND  BISTATIC  RCS  VALUES  READ  FROM  FILE 

C  RCSMON:  CALCULATED  MONOSTATIC  RCS  VALUES 

C  W:  FREQUENCY  IN  RADIANS /SEC 

C 

PROGRAM  RCSCONV 

DIMENSION  RCSE(721) ,RCSMON(721) ,A(15) 

C 

C  SPECIFY  VARIOUS  FORMATS  FOR  FILE  INPUTS/OUTPUTS 
100  F0RMAT(F6.2,1X,E13.6) 

C 
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OPEN(UNIT=50, FILE= ' rcsinon_v.plt' , STATUS= ' UNKNOWN ' ) 

DO  10  N=l,361 
NU= (N+19) /40+11 

OPEN (UNIT=ll,FILE='rcsm0_0. pit' , STATUS =' OLD ' ) 
OPEN (UNIT=12,FILE='rcsm0_10. pit' , STATUS =' OLD ' ) 
OPEN {UNIT=13 ,FILE='rcsm0_20. pit' , STATUS= ' OLD ' ) 
OPEN (UNIT=14,FILE='rcsm0_30. pit' , STATUS= ' OLD ' ) 
OPEN {UNIT=15,FILE='rcsm0_40. pit' , STATUS= ' OLD ' ) 
OPEN {UNIT=16,FILE='rcsin0_50. pit' , STATUS= ' OLD ' ) 
OPEN {UNIT=17,FILE='rcsm0_60. pit' , STATUS= ' OLD ' ) 
OPEN {UNIT=18 , FILE= ' rcsm0_70 .pit ' , STATUS= ' OLD ' ) 
OPEN {UNIT=19,FILE='rcsm0_80. pit' , STATUS= ' OLD ' ) 
OPEN (UNIT=20,FILE='rcsmO_90. pit' , STATUS= ' OLD ' ) 

READ(NU,*)  FILL,W 
DO  20  1=1,721 

READ(NU,*)  DEGREE , ROSE ( I ) 

CONTINUE 

IF(NU.EQ.ll)  THEN 
NB= (2* (N-1) +1) 

RCSMON(N) =RCSE(NB) 

ELSE  IF(NU.EQ.20)  THEN 
NB=-2* (361-N) +361 
RCSMON(N) =RCSE(NB) 

ELSE 

NB=2* (N- ( (NU-11) *40+1) ) + (NU-11 ) *40+1 
RCSMON{N) =RCSE{NB) 


END 

IF 

CLOSE 

(11) 

CLOSE 

(12) 

CLOSE 

(13) 

CLOSE 

(14) 

CLOSE 

(15) 

CLOSE 

(16) 

CLOSE 

(17) 

CLOSE 

(18) 

CLOSE 

(19) 

CLOSE 

(20) 

CONTINUE 

WRITE(50, 100)  FILL,W 
DO  50  1=1,361 

DEG=FLOAT(I-l) *0.25 
WRITE {50, 100)  DEG,RCSMON(I) 
CONTINUE 

DO  51  1=362,721 

DEG=FLOAT(I-l) *0 .25 
WRITE(50, 100)  DEG, RCSMON( 722-1) 
CONTINUE 


CLOSE  (50) 
STOP 


B.3  Additional  Codes 


The  following  code  listings  were  helpful  in  analyzing  and  visualizing  the  FVTD  results. 

B.3.1  Error  Calculations  and  Metrics 

The  following  MATLAB  code  calculates  the  difference,  in  dB,  and  the  mean-square- 
error  (MSE)  between  the  FVTD  and  MoM  results  and  the  cross-correlation  and  the  Fourier 
transform  of  the  results. 


%  This  program  uses  several  methods  to  calculate  various 
%  parameters  which  can  be  used  to  compare  FVTD  and  MoM  Results 
%  MATLAB  Code 
% 

%  The  various  methods  used  include: 

%  1.  Difference  in  dB  after  considering  phase  by  calculating 

%  the  minimum  error  between  the  FVTD  result  at  one  point  and 

%  several  surrounding  {+/-5  degrees)  MoM  results. 

%  2 .  Mean  Square  Error  (MSE)  using  FVTD  and  MoM  results  in 

%  square  meters.  Before  calculating  the  MSE,  phase  is  talcen 
%  into  consideration  in  the  same  way  as  for  method  1. 

%  3 .  Cross  Correlation  between  FVTD  and  MoM  results  in  square 

%  meters . 

%  4.  FFT  of  the  FVTD  and  MoM  results  in  square  meters. 

%  Load  the  data 
load  fvtd_hh.prn; 

FS=fvtd_hh; 

FS=FS(1;4:721)  ; 

FD=10.*logl0{FS) ; 

load  mom_hh.prn; 

MD=mom_hh ; 

MS=10 . ^ (MD. /lO) ; 

%  Method  1 

%  Calculate  difference  in  dB,  Adjust  Phase  (dB  vs.  theta) 

DDT= [0 : 1 : 180]  ; 
n=l; 

while  (n<=5 ) , 
m=l  ; 

A=100; 

while  (m<=6) , 

A=min(abs (FD(n) -MD(n- (m-6) ) ) ,A) ; 
m=m+ 1 ; 
end 

DDT(n) =A; 
n=n+l ; 
end 
n=177; 


%  FVTD  results  in  square  meters 
%  FVTD  results  in  dB 


%  MoM  results  in  dB 
%  MoM  results  in  square  meters 


142 


while  (n<=181) , 
m=6  ; 

A=100; 

while  (m<=ll) , 

A=min(abs  (FD(n)  -MD(n-  (in-6)  )  )  ,A)  ; 
in=m+l  ; 
end 

DDT(n) =A; 
n=n+l ; 
end 

n=6  ; 

while  {n<=176), 
m=l  ; 

A=100; 

while  (m<=ll) , 

A=min(abs  (FD(n)  -MD(n-  (in-6)  )  )  ,A)  ; 
m=m+ 1 ; 
end 

DDT(n) =A; 
n=n+l ; 
end 

DD=max ( DDT ) 

DDT=rot90 (DDT, 3) ; 

save  diffhh  DDT  -ascii  %  Save  difference  values  in  dB  to  a  file 
%  Method  2 

%  Calculate  MSE,  Adjust  Phase  if  Incorrect  (square  meters). ^2 

MSE=0; 

n=l; 

while  (n<=5) , 

MSE=MSE+(FS(n)-MS(n)  )  .''2; 
n=n+l ; 
end 
n=177; 

while  (n<=181) , 

MSE=MSE+(FS(n)  -MS(n)  )  .''2; 
n=n+l ; 
end 

n=6  ; 

while  (n<=176) , 
m=l  ; 

A=100; 

while  (m<=ll) , 

A=min(abs (FS (n) -MS (n- (m-6) ) ) ,A) ; 
m=m+ 1 ; 
end 

MSE=MSE+A.^2; 
n=n+l ; 
end 

MSE= (1. /181) . *MSE 
%  Method  3 

%  Compute  cross-correlation  (square  meters) 

XC=xcor (MS , FS , ' coef f ' ) ; 

%  Plot  the  cross-correlation 
theta= [0:1: (2*length(FS) -2) ] ; 
plot ( theta , XC ) 
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title {' Cross-Correlation  Between  FVTD  and  MoM  results,  HH') 

xlabel ( 'Theta  (Degrees) ' ) 

ylabel ( 'XCORR' ) 

axis{ [170  190  0.95  1] ) ; 

save  xcorrhh  XC  -ascii  %  Save  cross-correlation  values  to  a  file 
%  Method  4 

%  Calculate  FFT  of  FVTD  and  Mom  results  (square  meters) 

%  Remove  DC  value  first 
n=l; 

FST=0; 

while  (n<=181) , 

FST=FS (n) +FST; 
n=n+l ; 
end 

FST=FST. /181; 

FS=FS-FST; 

n=l; 

MST=0; 

while  (n<=181) , 

MST=MS (n) +MST; 
n=n+l ; 
end 

MST=MST. /181; 

MS=MS-MST; 

%  Use  Hamming  window 
FS=ham(181) . *FS; 

MS=ham(181) . *MS; 

%  Calculate  FFT 
FF=abs(fft(FS,1024)  )  ; 

MF=abs(fft(MS,1024)  )  ; 

%  Plot  the  fft's 
figure 

thetai=(l. 71024)  .*[0:1023] ; 
plot (thetai, FF) 
hold  on 

plot (thetai, MF,  ' -- '  ) 

title ('FFT  of  FVTD  vs.  MoM  Results,  HH') 

xlabel ( ' Cycles/Theta ' ) 

ylabel ( ' FFT ' ) 

axis( [0  0.1  0  0.06] ) ; 

save  fft_fhh  FF  -ascii  %  Save  FVTD  fft  values  to  a  file 
save  fft_mhh  MF  -ascii  %  Save  MoM  fft  values  to  a  file 


B.3.2  Scattered  Electric  Field  Movie 

The  listed  MATLAB  code  generates  a  movie  of  the  scattered  field  for  the  cone-sphere. 
The  plots  are  contour  plots  of  the  scattered  electric  field. 


%  Movie  of  Electromagnetic  Scattering  from  3-D  Object  (Cone-Sphere) 
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%  Movie  displays  contours  of  the  scattered  electric  field 
%  MATLAB  Code 

%  As  the  scattered  field  data  is  loaded  and  displayed  at  each  of  N  time 
%  steps,  snapshots  of  the  graphics  window  are  saved  in  a  large  matrix. 

%  The  total  storage  required  is  proportional  to  the  area  of  the 
%  graphics  window. 

%  Load  the  grid 
1=49; 

J=3  6  ; 

nframes=40 ; 

fid  =  f open ( ' X . grd ' ,  ' r ' )  ; 

[x, Count]  =  f scanf ( f id, ' %13e ' , [ I , J] ) ; 
fid  =  fopen{ 'y.grd' , 'r ' ) ; 

[y, Count]  =  fscanf ( f id, ' %13e ' , [I, J] ) ; 

%  Set  up  movie  frame 

fid  =  fopen ( 'plotSO . dat '  ,  ' r  '  )  ; 

[El, Count]  =  f scanf { f id, ' %lle ' , [ I , J] ) ; 

%  Plot  the  contour 
pcolor (x, -y, El) 
shading { ' interp ' ) 
colormap (jet) 

axis ( [-1.1  1.1  -1.1  1.7] )  ; 
axis ( ' square ' ) 
caxis ([0.0  0.3]); 

title (' Cone-Sphere  Scattered  Electric  Field  Contours  (0.869  GHz)') 
xlabel ( ' X  (m) ' ) 
ylabel ( ' y  (m) ' ) 
hold  on 

pcolor (-X, -y, El) ; 
shading ( ' interp ' ) 

M=moviein (nframes ) ; 

%  Record  the  Movie 
for  j=l: nframes 
if  j==l 

f id=fopen ( ' plotSO . dat '  ,  ' r ' )  ; 
end 

if  j==2 

f id=f open ( ' plotSl . dat ' , ' r ' ) ; 
end 

if  j==3 

f id=fopen ( ' plot52 . dat ' , ' r ' ) ; 
end 

if  j==4 

f id=f open ( ' plot53 . dat '  ,  ' r ' )  ; 
end 

if  j==5 

f id=f open ( ' plot54 . dat '  ,  ' r ' )  ; 
end 

if  j==6 

f id=f open ( ' plot55 . dat '  ,  ' r  '  )  ; 
end 

if  j==7 

fid= fopen ( ' plots 6 . dat '  ,  ' r  '  )  ; 
end 

if  j==8 

f id=f open ( ' plot57 . dat ' ,  ' r  '  )  ; 
end 

if  j==9 
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f id=fopen ( ' plot58 . dat ' , ' r ' ) ; 
end 

if  j==10 

f id=fopen ( ' plot59 . dat ' , ' r ' ) ; 
end 

if  j==ll 

f id=fopen ( ' plot60 . dat ' , ' r ' ) ; 
end 

if  j==12 

f id=fopen ( 'plot61 . dat ' , ' r ’ ) ; 
end 

if  j==13 

f id=f open ( ' plot62 . dat ' , ' r ' ) ; 
end 

if  j==14 

f id=fopen ( 'plot63 . dat ' , ' r ' ) ; 
end 

if  j==15 

f id=f open ( ' plot64 . dat ' , ' r ' ) ; 
end 

if  j==16 

f id=f open ( ' plot65 . dat ' , ' r ' ) ; 
end 

if  j==17 

f id=fopen ( ' plot66 . dat ' , ' r ' ) ; 
end 

if  j==18 

f id=f open ( ' plot67 . dat  ’  ,  '  r ' )  ; 
end 

if  j==19 

f id=fopen { ' plot68 . dat '  ,  ’  r ' )  ; 
end 

if  j==20 

f id=fopen ( ' plot69 . dat ' , ' r ' ) ; 
end 

if  j==21 

f id=f open { ' plot7 0 . dat ' ,  ' r ' )  ; 
end 

if  j==22 

f id=fopen ( ' plot71 . dat '  ,  ' r ' )  ; 
end 

if  j==23 

f id=fopen { ' plot72 . dat ' , ' r ' )  ; 
end 

if  j==24 

f id=fopen ( ' plot73 . dat '  ,  ' r ' )  ; 
end 

if  j==25 

f id=fopen ( ' plot74 . dat ' , ' r ' ) ; 
end 

if  j==26 

f id=fopen ( ' plot7  5 . dat '  ,  ’ r ' )  ; 
end 

if  j==27 

fid=fopen{ 'plot76 .dat ' , ' r ’ ) ; 
end 

if  j==28 

f id=f open ( ' plot77 . dat ’ , ' r ' ) ; 
end 

if  j==29 
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f id=f open { ' plotV  8 . dat ' ,  ' r ' ) ; 
end 

if  j==30 

f id=fopen ( ' plot? 9 . dat ' , ' r ' ) ; 
end 

if  j==31 

f id=fopen ( ' plotSO . dat ' , ’ r ' ) ; 
end 

if  j==32 

f id=fopen { 'plotSl . dat '  ,  '  r '  )  ; 
end 

if  j==33 

f id=f open ( ' plot82 . dat ' , ' r ' ) ; 
end 

if  j==34 

f id=f open ( ' plot83 . dat '  ,  '  r  '  )  ; 
end 

if  j==35 

f id=f open ( ' plot84 . dat '  ,  '  r  '  )  ; 
end 

if  j==36 

f id=f open { ' plot85 . dat '  ,  '  r '  )  ; 
end 

if  j==37 

f id=f open ( ' plot86 . dat '  ,  ' r ' )  ; 
end 

if  j==38 

fid=fopen{ 'plot87 .dat' , 'r ' )  ; 
end 

if  j==39 

f id=fopen { 'plot88 . dat '  ,  ' r ' )  ; 
end 

if  j==40 

f id=fopen { ' plot89 . dat '  ,  '  r  ’  )  ; 
end 

[El, Count]  =  fscanf {fid, ' %lle' , [I, J] ) ; 

%  Plot  the  contour 

pcolor (x, -y, El ) 

shading ( ' interp ' ) 

axis( [-1.1  1.1  -1.1  1.7] ) ; 

axis ( ' square ' ) 

caxis ([0.0  0.3]); 

title (' Cone-Sphere  Scattered  Electric  Field  Contours  (0.869  GHz)  ') 

xlabel { ' X  (m) ' ) 

ylabel ( 'y  (m) ' ) 

pcolor  (-X, -y,  El)  ,- 

shading ( ' interp ' ) 

M{ : ,  j ) =getframe; 
end 

%  Play  Movie 
movie (M, 5,4) 


B.4  FVTD  FORTRAN  Code 

The  entire  code  can  be  obtained  from  Dr.  J.S.  Shang  of  Wright  Laboratory  or  the  author. 
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